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ABSTRACT 

We derive new equations using the mixed-frame approach for one- and two-dimensional (axisymmet- 
ric) time-dependent radiation transport and the associated couplings with matter. Our formulation is 
multi-group and multi-angle and includes anisotropic scattering, frequency(energy)-dependent scat- 
tering and absorption, complete velocity dependence to order v/c, rotation, and energy redistribution 
due to inelastic scattering. Hence, the "2D" realization is actually "6 1/2" -dimensional. The effects 
of radiation viscosity are automatically incorporated. Moreover, we develop Accelerated-Lambda- 
. Iteration, Krylov subspace (GMRES), Discontinuous-Finite-Element, and Feautrier numerical meth- 

ods for solving the equations and present the results of one- dimensional numerical tests of the new 
formalism. The virtues of the mixed-frame approach include simple velocity dependence with no 
.~ . velocity derivatives, straight characteristics, simple physical interpretation, and clear generalization 

to higher dimensions. Our treatment can be used for both photon and neutrino transport, but we 
\^ \ focus on neutrino transport and applications to core-collapse supernova theory in the discussions and 

examples. 

Subject headings: multi-dimensional radiation transport, radiation hydrodynamics, numerical meth- 
ods, supernovae, neutrino transport 
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Many phenomena in the Universe must be addressed using the tools of radiation transport and radiation hydro- 
dynamics to achieve a theoretical understanding of their character. Supernova explosions, gamma-ray bursts, star 
formation, planet formation, nova explosions. X-ray bursts, Luminous-Blue- Variable (LBV) outbursts, and stellar 
^ ' winds are all time-dependent fluid flow problems in which radiation plays a pivotal, often driving, role. In most 
circumstances, the radiation is photons, but in studies of the core-collapse supernova mechanism the radiation is neu- 
Q ,. trinos. In either case, time-dependent techniques to address radiation transport and the coupling of radiation with 
$^ matter are of central concern to the theorist whose goal is explaining the transient, dynamical phenomena of the 
^ \ Cosmos. However, spherically-symmetric algorithms for radiation transport and radiation hydrodynamics, generally 
^ I necessary to achieve a flrst-order understanding, are often not sufflcient and multi-dimensional approaches are called 
for. These are not easy, not only to formalate, but to implement. Nevertheless, multi-D radiation transport is emerging 
as a necessary tool in the theorist's toolbox and computers to evolve the associated equations are becoming available 
to a wider cohort of researchers. 
^ ■ In this paper, we derive new equations using the mixed-frame approach for one- and two-dimensional radiation 
5^ , transport and the associated coupling with matter, significantly extending the pioneering work of Mihalas & Klein 
(1982) to include anisotropic scattering, frequency(energy)-dependent scattering and absorption, complete velocity de- 
pendence to order ti/c, rotation, and energy redistribution due to inelastic scattering. Moreover, we develop algorithms 
for solving these equations and present the results of one-dimensional numerical tests. Hence, we provide both the 
new formulation and appropriate numerical techniques to solve it. The mixed-frame approach, in which the radiation 
quantities are defined in the laboratory (Eulerian) frame and the matter and coupling quantities are defined in the 
comoving frame, has largely been neglected by the radiation transport and atmospheres communities because of their 
focus on line transfer. The desire to include narrow spectral lines and to handle Doppler shifts into and out of those 
lines necessitates many spectral bins and a huge number of angular bins to ensure the lines are resolved on the compu- 
tational grid. As a result, most dynamic atmosphere and radiation studies are done using the comoving (Lagrangian) 
equations of radiation transport. For instance, the core-collapse supernova community, which has been at the forefront 
of radiation hydrodynamic developments in astrophysics, has inherited this formulation for their treatment of neutrino 
transport. 

However, many radiation-hydrodynamic problems do not require an exquisite treatment of spectral line transport, 
but a good treatment of continuum transport. The core-collapse supernova problem is one such case. The monochro- 
matic opacities and emissivities of neutrinos are overwhelmingly smooth functions of neutrino energy. Given this, the 
mixed-frame formulation is ideally suited for supernova theory. Its virtues vis a vis the comoving-frame approach to 
the solution of the Boltzmann equation and its related moment equations are numerous: 1) Even in one-dimension, 
instead of requiring ~20 velocity-dependent terms (Buras et al. 2006) on the left-hand(streaming) side of the Boltz- 
mann/transport equation, many of which involve spatial velocity derivatives, there are no such terms on the left-hand- 
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side in the mixed frame approach and only one grouped linear term (to 0{v/c)) on the right-hand(source) side; 2) 
There are no terms with derivatives of the velocity. Therefore, the characteristics of the associated transport equation 
are all straight lines. Furthermore, there is no need for the monotonicity in the velocity field required by some implicit 
solvers; and 3) The mixed-frame method is easily generalized to two and three dimensions, and the associated solvers 
are straightforward (though more expensive) extensions of those employed in ID. Note that much has been made of the 
importance of velocity-dependent terms in the transport equations for the calculation of the neutrino energy deposited 
in the "net-gain" (Wilson 1985) region. We show that the mixed-frame approach provides the most straightforward 
perspective from which to understand the physics of this effect and that even its sign is frame-dependent. In particular, 
in the mixed-frame formulation the velocity-dependent term augments the net gain during the stalled-shock phase, 
while in the comoving frame formulation the corresponding terms reduce it. Hence, any statement concerning the 
importance of such terms is very frame- and treatment-dependent. 

Though the mixed-frame equations depend simply upon velocity, the Lorentz transformations that are the core 
of this formulation introduce frequency (energy) derivatives of the radiation moments. These couple adjacent energy 
groups and might have compromised implicit algorithms that parallelize in group (Livne et al. 2004; Walder et al. 
2005; Burrows et al. 2006; Dessart et al. 2006ab; Ott et al. 2006ab). However, we show that such terms can be handled 
semi-implicitly, with the logarithmic derivatives of the moments with respect to energy handled explicitly during an 
otherwise implicit solve. In this way, processors performing updates on only a single group are not coupled during 
the iteration to adversely affect parallelization and scalability. Furthermore, a similar semi-implicit tactic works for 
inelastic scattering terms, since these are sub-dominant in the context of core collapse (Thompson, Burrows, & Pinto 
2003). Note, however, that parallelization in energy groups is viable for 2D calculations, but not for 3D. An entire 2D 
hydro grid can now reside on a single processor, but an entire 3D hydro grid requires spatial domain decomposition onto 
many processors and parallelization in only encirgy groups is not yet viable. Furthermore, given the seven-dimensional 
nature of 3D radiation transport (3 space + 2 angles + 1 frequency /energy -|- 1 time), 3D is not yet computationally 
feasible for astrophysical simulations. Therefore, we focus in this paper on ID and 2D mixed-frame formulations. 

However, a 2D, mixed-frame, azimuthally-symmetric formalism is still six-dimensional (lacking one spatial dimen- 
sion and a hemisphere) and is particularly straightforward and powerful. Our approach includes rotation (and, hence, 
is actually "6 1/2" -dimensional). The Eddington factor becomes an "Eddington tensor" with five independent en- 
tries. The zeroth- and first-moment equations are accurate approximations to the full equations that are closed with 
Eddington factors that can be calculated at each timestep to provide the full solution, or every N steps to provide 
an excellent, fast, though approximate, solution. Since velocity dependence is included in the multi-D context, the 
effects of radiation viscosity (neutrino or photon) are automatically incorporated into any scheme that includes the 
hydrodynamics with the radiation force and energy couplings. 

Though our treatment can be used for both photon and neutrino transport, it was devised with neutrinos in mind. 
As a result, much of the discussion and all the tests we perform are in that context. Therefore, a short review of 
the transport schemes employed to date in supernova modeling is in order. In the realm of ID (spherical) neutrino 
transport. Bowers fc Wilson (1982), Bruenn (1985), Wilson (1985), and Mayle, Wilson, & Schramm (1987) employed 
multi-group diffusion codes with flux limiters and did not address angle-dependent transport. Later, various groups 
achieved a multi-group/multi-angle neutrino Boltzmann capability (Mezzacappa & Bruenn 1993; Rampp & Janka 
2000; Mezzacappa et al. 2001; Thompson, Burrows, & Pinto 2003; Liebendorfer et al. 2001ab, 2004, 2005). In 2D 
simulations, LeBlanc & Wilson (1970) constructed a gray, flux-limited diffusion code, that nevertheless calculated a 
two-component vector flux. Herant et al. (1994) fielded a gray diffusion code with a simple matching to the free- 
streaming regime. Burrows, Hayes, & Fryxell (1995) developed a gray diffusion code that calculated different solutions 
for the different 6 angles (in spherical coordinates) in a 2D hydro simulation, but calculated at each angle a spherical 
model that employed the matter profiles along that one ray. They did not calculate lateral transport in the angular 
direction. This is the so-called "ray-by-ray" approach that is still used by some groups today. Rampp & Janka 
(2000,2002), Janka et al. (2005ab), and Buras et al. (2006) have updated the ray-by-ray method using sophisticated 
ID Boltzmann transport along each ray, but no lateral transport (though they do include in the hydro the effects 
of lateral radiation pressure and lepton number transport). The transport is calculated in the comoving frame and 
mapped to the Eulerian frame of their PPM hydrodynamics. This prescription is viable only if the flow maintains 
rough sphericity and the core does not move off the center of the grid (as when kicks are imparted to the protoneutron 
star). However, the ray-by-ray approach is not real 2D transport. The ZEUS-2D radhydro code of Stone, Mihalas, and 
Norman (1992), and its update by Hayes and Norman (2003), solve the zeroth- and flrst-moment equations of the gray 
transport equation and use a short-characteristics solution of the static transfer equation to obtain the second moment 
for closure. These codes arc realizations of the variable Eddington factor method and avoid some of the pitfalls of flux 
limiters, but since they are gray and not multi-group, they are of limited utility for modern supernova simulations. 

All the above multi-group codes were formulated in the comoving frame and none of the 2D variants was simulta- 
neously 2D, multi-group, and time-dependent. A 2D, multi-group, flux-limited capability has recently been developed 
by Swesty & Myra (2005ab, 2006), but they treat the inner core in ID, not 2D, and to date have not pubhshed long- 
duration core-collapse and post-boimcc simulations. Cardall, Lentz, & Mezzacappa (2005) and Cardall & Mezzacappa 
(2003) have derived the general form of the neutrino transport equations including higher-order terms in v/c and 
general relativity, but have not yet developed working implementations. The first bonafide 2D multi-group, multi- 
angle, time-dependent capability was achieved by Livne et al. (2004) using the implicit, Arbitrary-Lagrangian-Eulerian 
(ALE) code VULCAN/2D (with remap), for which the radiation field was defined in the laboratory frame, but this 
code is not fast, does not include the Doppler shifts due to the velocity field (though it does include advection), and 
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does not include energy redistribution. However, its multi-group, flux-limited diffusion variant is 2D in the entire 
computational domain and much faster than its multi-angle version. With it, a number of multi-group, fully-2D, 
radiation/hydrodynamic investigations have been possible, with and without rotation (Walder et al. 2005; Burrows et 
al. 2006; Ott et al. 2006ab; Dessart et al. 2006ab). However, as flops become cheaper and computer speeds improve, 
one will need to do better. This is what motivates the present paper and the development of the new 2D mixed-frame, 
implicit, multi-group, multi-angle algorithm and its two-moment variants. In the context of hydrodynamics, this code 
has been christened BETHE^. 

In ii2.1l we present the transport equation and its general formulation, but quickly in t |2.2l derive. using the appropri- 
ate Lorentz transformations, the equations of radiative transfer in the mixed-frame. This section contains our central 
analytic results. In W2.'d\ we follow with derivations of the associated moment equations. Then, in ti2.4l we digress into a 
discussion of the frame transformations of the source terms on the right-hand-sides of the transport equation and of the 
associated matter energy and momentum equations. Care in these matters is important to ensure global energy and 
momentum conservation to 0{v/c) and that one employs the correct radiation source terms in the hydro equations. 
Different realizations of the hydro equations are possible and the form of the radiation energy source term depends 
upon the form of the hydro energy equation. For instance, the matter energy equation can be in first-law format 
(purely Lagrangian with comoving-frame time derivatives) or can include the kinetic energy explicitly (with Eulerian 
partial time derivatives). For these two formulations, the radiation source terms are different. In ^ we present the 
mixed- frame formalism in cylindrical (axisymmetric) coordinates and in 21 present the mixed- frame formalism in 
spherical coordinates. Then, in [j3 derive the mixed-frame equations in spherical symmetry and explore various 
solution techniques, including Accelerated-Lambda-Iteration (ALI), a tridiagonal approximate operator, the use of the 
Krylov subspace algorithm GMRES, the Discontinuous-Finite-Element (DEE) method, and the Feautrier scheme. In 
H5.5I we derive the associated moment equations, provide the matrix representation, and introduce sphericity factors. 
In ijni we derive and discuss procedures for the implicit coupling of matter with radiation, including an ALI treatment 
and the linearization of the energy and composition (Ye) evolution equations. We follow this in fJ7|with a series of 
numerical tests using the spherical formulation of the mixed- frame transport equations. These include stationary solu- 
tions f ii7.1(l and a time-dependent cooling calculation of an idealized protoneutron star with implicit radiation-matter 
coupling fi i7.2(l . Section |H1 provides a summary and the Appendix contains a derivation of the mixed- frame treatment 
of inelastic energy redistribution. 

2. FORMULATION 
2.1. Transport Equation 

The Boltzmann transport equation for the neutrino occupation probability / is (Mihalas & Mihalas 1984): 

i ^ + (n • V) / = Cth [/] + Ces [/] + Cnes [/] , (1) 
c ot 

where / is the neutrino occupation probability, n the unit vector in the direction of neutrino propagation, Cth is 
the collision integral (net source term) for "thermal" creation and destruction of neutrinos (emission and absorption) , 
mostly due to charged-current processes; Ces is the coUisional integral for elastic scattering of neutrinos; and Cnes is 
the collision integral for inelastic scattering (such as neutrino-electron scattering). We rewrite the transport equation 
using the specific intensity, /, which can be written in terms of / as: 

where v is the neutrino energy, h is Planck's constant, and c is the speed of light. 

It is customary to formulate the combined Cth + Ces contributions in terms of absorption and emission coefficients. 
The transfer equation then reads 

i ^ -f n • /(i., n) = (^^ n) + r;- (z., n) - [a^ (i^, n) + a(z., n)] /(z., n) , (3) 

where k is the true absorption coefficient, a is the scattering coefficient, is the thermal emission coefficient, and 
T]^'^ is the scattering part of the emission coefficient. This equation has the same form in all frames. 

Inelastic scattering is somewhat complicated, but fortunately it is usually small compared with the thermal and 
elastic scattering source terms. We postpone a detailed investigation of inelastic scattering to a future paper, but 
provide its formulation in the mixed-frame formalism in the Appendix. 



2.2. Mixed-Frame Formulation 

In the mixed-frame approach, the material properties (absorption and emission coefficients) of the right-hand side 
of eq. (O are expressed using the comoving frame, while the specific intensity and the left-hand side are expressed 
in the inertial, observer's frame. This approach was first suggested in the context of photon transport by Mihalas 
& Klein (1982, hereafter referred to as MK). In this paper, we will generalize their approach by allowing for energy- 
dependent anisotropic scattering, as well as for non-coherent, inelastic scattering (Appendix). While the mixed- frame 
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formalism has only a limited applicability for photon line transport due to a large variation of spectral line opacity 
as a function of photon energy, the mixed- frame approach is very well suited to neutrino transport, in which all the 
relevant neutrino-matter interaction cross-sections are smooth functions of neutrino energy. 

Denoting by subscript quantities in the comoving frame, the Lorentz transforms of the photon/neutrino energy 
and direction are 



1^0 = ^1 



and 



no = {v/vo) 
To 0(v/c), we have the following expressions 



V 

n I 7 



7 n • V 

7+1 c 



V 

c 



(4) 
(5) 

(6) 
(7) 

(8) 
(9) 

(10) 

Both absorption coefficients in the incrtial frame are expressed through the comoving frame coefficient and its 
derivative, exactly as in MK 

n • V 

K{i^,n) Ko(i^) 



/ V \ / n • V 

no = (i^/i^o) [n- -j =n[l^ — 

and the absorption (k) and scattering (cr) coefficients transform as 
and 

(^{v^^) {i'o/v)(^o{vq) . 
The emission coefficient in eq. (jSJ transforms as 



(11) 



which follows from eqs. © and © and a Taylor expansion k.o{vo) = '*o('^) + {1^0 ^ j^)9ko/9i^. The transformation of 
a is analogous: 



The thermal emission coefficient, as given by MK, is: 



th/ 



n • V 



dag 
dv 



du 



where we assume that the thermal emission coefficient is isotropic in the comoving frame. 
We assume that the comoving-frame elastic scattering emission term is given by 

r/g'=(i/o,no) = j> duj'ah[vo,-n'o) go{-n'Q,no) , 



(12) 



(13) 



(14) 



where the primed quantities refer to the properties of the absorbed photon/neutrino and is the scattering phase 
function in the comoving frame. 
In the following, we assume that the scattering phase function is in a simple form. 



5o(no,no) = 1 + (5 n^ • no 



(15) 



For elastic neutrino-matter scattering this is an excellent approximation. The scattering emission coefficient transforms 
according to eqs. (fTH)! and (fl^ as 



ry-(i.,n)= - 



V \ cro(t'o) 



47r 



where X is the integral term of eq. 1)14(1 . The specific intensity transforms as 



/o(^^,n;,) = (^)'/(^',n'). 



and the element of the solid angle as 



(16) 

(17) 
(18) 
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We express the specific intensity at v' through the Taylor-series expansion around 

I{v\ n' = n' + I{v, n') + v — [ , 19 

av ov \ - 

and the cosine of the scattering angle, again to 0{v/c), as 

/ , / / , / n' ■ V n • , 1 , , 

nf, • no = n' ■ n + (n' • n - 1) + . (20) 



A final step is the connection between the conioving and laboratory frames, which consists in accounting for the Taylor 
expansion of (Tq to transform the energy from vq to i>: 

/ \ / \ n • V dao , , 

croWo) = cro(iy) (21) 

c oiy 

To avoid confusion, we stress that we do not need here the transformation equation (|12|) that describes the transfor- 
mation from the inertial-frame a to the comoving- frame cto- Here, we already have the comoving- frame ao, and all we 
need is to transform the energy, which is exactly what eq. H21(l expresses. 
After some algebra, we obtain for the transport equation in the mixed-frame formalism: 

Id \ 

-— + n-Vj I{iy,n) = roo(i^, n) + roi(i^, n) + rio{iy,n) + rii(i/, n) , (22) 

where the subscripts refer to terms which are 0-th or 1-st order in v/c and S, respectively. The terms of the r.h.s of 
eq. (I22|l are given by 

ron=Vo'-{i^Q + <='o)I{'',ri)+(ToJ, (23) 

roi=<7o6Wn,, (24) 



+njW-' 



9 In (To 



and 



rii = CTo (5 



9 In J/ 

rijW^ J + UjWk 



rio = r 
I{v,n) + aoJ ( 2 

aini^J'^l 



d in 77q 



th 



9 In (To 9 In J 
9 in d In v 



Ko 1 



9 In Ko 



d\ni' 
aoWjH^ I 1 



dlnv 



aoSW 



n) 

d\nW 
dhiv 

9 In (To d In 



d\a.v 



91n ^' 



(25) 



(26) 



where, again, we omit explicit indication of the energy dependence of most quantities. Here we use the convention that 
one sums over repeated indices, although our use of subscripts or superscripts is arbitrary and does not necessarily 
denote covariant or contravariant components. To simplify the notation, we have introduced the normalized velocity. 



V 

w = — 

c 



Here, the usual moments of the specific intensity are defined by 

<j> I duj /47r , 



47r 



H 



— = (f Induj/A-K . 



cP 

47r 



<j) I nn du /4,iT , 



(27) 

(28) 
(29) 
(30) 



where E, F, and P are the radiation energy density, flux, and pressure tensor, respectively. 

2.3. Radiation Moment Equations 
Integrating transfer cq. Q over angles with the source/sink terms given by H23(l - H26|l . we obtain the 0th- and 



1st- moment equations (written here in Cartesian coordinates) 

1 dJ dW 
c dt dx^ 

and 



(31) 



1 dK^i 



c dt 



= -(ko + o•tr)i^' + w'f^o + y CTq J [2-5- 
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3 91ni 



(32) 
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In eqs. (I31|l - (|33|) . we have introduced the so-called transport cross-section, 
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CTtr = (To 1 - 



and have set 



and 



Kq = Kq [l 



d In Kq 
1 



Vo 



Co = ctq I 1 

9 In j]^ ' 
9 In J/ 



9 In (To 
91ni^ 



To close the system of moment equations, we introduce the Eddington tensor: 

m = 1^^ 

^ " J ' 

so that the Ist-moment equation is written 

IdH' d 



where 



Wi 



(To 2 - ,5 



c dt dxi 

d In (To 9 In J 
9 In 9 In 



Ko + (To 1 



d In i 



3 aini 



3 aini 



(33) 

(34) 
(35) 
(36) 

(37) 

(38) 
(39) 



2.4. Hydrodynamical Equations and Energy Coupling 

We may write moment equations H31|l and (|32|l as equations for the energy and momentum of the radiation field, 
namely 

'^'^-cG\ (40) 



and 



dt dx 
1 OF' dP'^ 



c2 dt 



(41) 



where here E, F, are P are the energy-integrated radiation energy density, flux, and stress tensor, respectively, and 
G° and are the components of the four- force density vector (Mihalas & Mihalas 1984), in the inertial frame, which 
are given in our mixed-frame formalism by the trivial modification of the r.h.s. of eqs. (31) and (32), viz. 



/"OO 

cG° = 47r / [ko J - vl^ - W] dv , 

"^O 



and 



cG' = 47r 



[(ko + CTtr) H' ~ w'lio - e J] dv . 



(42) 



(43) 



The equation for overall energy conservation of the radiating fluid (i.e., containing both the matter and neutrino 
energy), correct to 0{v/c), is given by (MK; Mihalas & Mihalas 1984): 



^1 ^1 



(44) 



where e is the specific internal energy of the fluid, p is the fluid pressure, and /' is the external force on the fluid. This 
equation can also be written as an total energy for the radiating flow. 



p2.ie + vy2) + ^{pv^)^rHr- 



dE OF' 
'dt^'d^ 



= zhP + cG° , 



(45) 



and analogously the momentum equation 
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To obtain the gas-energy equation, one first writes an equation for mechanical energy which is obtained by multiplying 
eq. (|46|l by Vi, and subtracts it from the total energy equation H45|l . The equation for mechanical energy, to 0{v/c), 

Divy2) „ dp ( I dF' dP'^\ 



The resulting comoving-frame gas-energy equation reads 

De , Dil/p) 



Dt ' ^' Dt 



= cG" - ViG' . (48) 



This is the appropriate equation for updating the fluid temperature and is a statement of the flrst law of thermody- 
namics. If we were to leave density fixed, we could write the energy equation as an equation for temperature, 

DT 

pGv^^cG" ~v,G\ (49) 

where T is the temperature, and Gy the specific heat at constant volume. 

Equation 149|) contains the components of the inertial-frame four-force density vector, which are given by eqs. H42|l 
and H43|l that in turn were derived by using the mixed-frame formalism. However, cG^ and cG* contain terms that are 
proportional to the scattering coefficient ao- Since the elastic scattering should not contribute to the energy balance 
in either inertial or comoving frames, one should make sure that such terms cancel exactly. This is very important in 
the context of neutrino transport in supernovae, since in a low-temperature, low-density regions beyond the shock the 
scattering coefficient ao may be larger by many orders of magnitude than the true absorption coefficient kq, and one 
can, thus, easily introduce spurious terms in the energy balance because of rounding errors. 

It is, therefore, instructive to derive the right-hand-side of the energy equation in a different way, in which we eliminate 
the scattering terms analytically. The idea is first to write down the four-force density vector in the comoving frame, 
where it is easy to formulate, and then to perform a Lorentz transformation back to the inertial frame. 

The four-force density vector in the co-moving frame is simply written as 

cGg==^ ^[xo(i'o)^(i'o,no) -7?o'(f^o) -'7o''(f'o,no)]7io di^odwo, (50) 
where the scattering term in the comoving frame is simply given by 

'7o''('^o, no) = <Jo{iya) / f Io{vo, 5o(no, no) {(Ljq/At:) 



= j> no)(l + (5no • no) {(Iuj'q/Att) 

= (To(i'o) [Jo(j^o) + '5Ho(i^o) • no] . (51) 

Therefore, 

/•oo 

cGl = A^ [«o(i^o)Jo(j^o)-^*('^o)]di'o, (52) 



and 

cGl = 47r / [ko(j/o) + (1 - <5/3)ao(i^o)] Hl{vo) dvo = An I [kqM + (^tM] HqM dv^ . (53) 



We can now easily transform the components of the four-force density to the inertial frame using the standard Lorentz 
transform [c.f., Mihalas & Mihalas 1984; eqs. (91.22)], which to 0{v/c) become: 

cG" ^ cG^o + ^^Go , (54) 

and 

cG' = cGl + v'Gl . (55) 

The inertial-frame four-force density vector is thus (dropping explicit indication of its dependence on the comoving- 
frame frequency vo) 

/•oo 

cG° = 47r / [kqJq - rf^ + w,{no + CTtr)i?^(j^o)] dvQ , (56) 
Jo 

and 

/•oo 

cG^ = Att / [(^0 + <ytM + w'M - Vo'] dvo , (57) 



and, thus, the right-hand-side of the material internal energy equation reads 

nOO 

cG° - ViG' = cG° - w, cG' = 47r / [koJqM - Vo'] dv^ -f 0{v'^/c^) . (58) 
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We see that the second term of eq. (I56|) and the first term of eq. (|57|l exactly cancel, as can be expected on physical 
grounds, and eq. (|5^ is the result. 

As a final step, we now have to transform the comoving frame moments of the specific intensity back to the inertial 
frame (while leaving the material quantities in the comoving frame untouched in keeping with the general philosophy 
of the mixed- frame approach) . 

Specifically, 

/"OO 



dv 



dhw 

hq{u)J{v) — a{v) ■ w ( 2ko 



kq{vq)Iq{i'q, no) dvf) duo 
(1 — 2n • w) n) dv dto 



dlniy 



dv , 



(59) 



and, thus, the appropriate right-hand-side of the energy equation in the inertial frame, written in the mixed-frame 
formalism (in which the material quantities are in the comoving frame, while the radiation quantities and the energy 
are in the inertial frame) is: 



cG° - ViG' = An 



91n^ 



dv . 



(60) 



We now have to show that we obtain the same expression if we use cG° and cC given by eqs. and (|^ . Using 
eqs. (1^ and we have 



1^0 + CTtr 



aino-Q SlniJ* 



and 



cG' = 47r 



9 In dlnv 
Ko + (Ttr)H' dv +0{w'), 



dv , 



(61) 
(62) 



and, thus, to 0{v/c), we obtain 

cG° - cG' = 47r / {hqJ - rjl"") dv - in I w^H' ( 2ko 



-47r 



ftr 



a In (To dlnW 



dlnv 9 In 



dlnv 



dv 
dv . 



(63) 



The last integral in eq. H63|) is equal to zero, as can be easily shown by integrating its first two terms by parts. 
Equation is thus, indeed, identical to eq. (|^ . Here we assume, consistently with the previous formalism, that 
the incoherence parameter S is independent of energy. One could easily develop a formalism in which one can account 
for the energy dependence of S (containing additional terms proportional to w^6 dlnS/dlnv), but we consider such a 
complication unnecessary. 
Next, we consider the electron fraction equation. In the comoving frame, it is given by 



pNa- 



DY,. 



Dt 



roc 

•Jq 



dvp 

1^0 



(64) 



where the sum extends over the neutrino species, and Si = —1 for Vg neutrinos, Si = 1 for v^ neutrinos, and Si = 
for all other neutrino species; Ny^ is the Avogadro's number. The right-hand-side of eq. (|64|l is easily expressed in the 
mixed frame using an analogous procedure as that used in eq. I59|l . which differs only by the occurrence of the term 
dvo/vQ instead of dvQ. The resulting equation is: 



pNa- 



Dt 



-An^s, 



KqJ 



d In h 



dv 

V 



Next, we consider the momentum equation. It was already given by eq. H46|l . namely 

Dv' _ . dp 



P 



-I- G' - G° , 



(65) 



(66) 



The radiation-interaction term, which can be viewed in Cartesian coordinates as a gradient of radiation pressure, 
—dpll^^/dx^ = G* — w^G'^, is written in eq. H66|l in the inertial, Eulerian, frame. To transform it to the mixed frame, 
we use the same strategy as before to obtain the net radiation heating term in the energy balance equation: we can 
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either use expressions for and G" given by eqs. H42|l and H43() . which are already expressed in the mixed- frame 
formalism and perform necessary integrations over energies, or we can use a simple expression for the comoving-frame 
four-force vector, eq. H53|) . and express the comoving-frame moment Hq(vq) through the inertial-frame moments. 
Obviously, both approaches have to yield the same result, which is: 



An 



{kq + (Ttr) H' - w' J + W. 



d{PJ) 



dv . 



(67) 



Finally, we stress the following feature of our formalism. The material equation for the conservation of the total energy 
and momentum, as well as the electron fraction equation, were written in the Eulerian frame, where the radiation- 
interaction is also in the Eulerian frame with, however, interaction coefficients (i.e., the absorption, emission, and 
scattering coefficient) in the comoving (Lagrangian) frame. If the overall hydro scheme is Eulerian, our present scheme 
is obviously consistent with it and we would use eq. (|42|) . If the overall hydro scheme is Lagrangian, the radiation- 
interaction terms are generally different. For instance, the right-hand-side of the gas-energy equation expressing the 
Ist-law of thermodynamics is given by cG° — WiG^ in Eulerian frame quantities [see eqs. H48(l and H60|l ]. while it is 
given by cGq, eq. H52|l . in comoving-frame quantities, which is formally different. 

3. MIXED-FRAME RADIATION EQUATIONS IN CYLINDRICAL GEOMETRY 

Here, we assume cylindrical geometry with azimuthal symmetry. The coordinates are r, z, and (j). We assume that 
the rest- frame material properties (temperature, density, opacity, emissivity, etc.) depend only on r and z, while the 
velocity field has a non-zero ^-component. In other words, we allow for rotation. Such an approach is sometimes called 
the "2 1/2-D" case in hydrodynamics. 

The unit vector in direction n is 



n = sin 6 cos + cos 9 + sin 9 sin ij^ ga 



(68) 



Here 9 is the polar angle measured from the positive z-direction and ^ is the local azimuthal angle, such that i]j = Q 
is in the local positive r-direction. 
In component form, the first and second moments are given by 



Hr = <j> sin 9 cos^ n)(i(jj/(47r) . 



cos9 n)duj/{4Tr) , 



and 



= <j) sinf? BUYip I^u^ n)(ia;/(47r) , 
Krr = <j> sin^ 9 cos^ -0 n)dcj/(47r) , 



cos 9 sinfl cos ip n)dLu/{4:Tr) , 
Kr^ — K^r — 'j> sin'^ 9 sinV' cosij} n)(iijj/(47r) , 

K^^ = <j> cos^ 9 1{v, n)da;/(47r) , 
K^z — Kztj, — <j) cos9 sin0 sm.^ I {y,i\)dw / {^Att) , 



and 



K. 



00 



sin^ 9 siv? ij} I{v,n)duj/{'i'K) , 



which, due to symmetry and the trace condition 

Krr + Kzz + Kss ~ J , 



(69) 
(70) 
(71) 

(72) 
(73) 
(74) 
(75) 
(76) 

(77) 
(78) 



leaves five, instead of nine, independent, non-zero components of the tensor K. 

The radiative transfer equation, in cylindrical coordinates and with azimuthal symmetry, written in the conservative 
form, now becomes 



1 dl^ „ dljj sin 9 cos ijj d , ^ . sin 6* 5 , . , ^ , 

— ^ +C0S6' — H Try^^'^l ^(sinV'4) 

c at oz r or r oip 

= roo + »'oi + ^10 + ni , 



(79) 
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where the right-hand side is given by eqs. (|23|l - ()26fl written in component form, viz. 

roo = vl^ - {ko + c^o) + croJ , (80) 
roi — ctqS {Hr sin 9 cos tp + Hz cos 9 + sin 9 sini/;) , (81) 

rijW-' = sin 9 cos ij^Wj. + cos 9wz+ sin sin , (82) 

and 

TLjWkK-'^ = sin0 cosi/i {wrK„ + WzK^z + w^Kr^) 

+ C089 {WrKrz + W^Kzz + W^Kz4>) + S\Ii9 fi\n%}) {WrKr 4, + WzKz^ + W^K^^). (83) 

The corresponding moment equations read: 

1 1 c^i? 

-7^ + -Tr(^^fr.) + ^ = ??o^ - KoJ + S,F, + ^zHz + 50i^0 , (84) 
c at r or oz 

2^ \ Q Q \ j j 

-^rr + -Tr(r/rr J) + Tr(/rz^) — —J= - {kq + atr)Hr + Wrrjo + ^rJ , (85) 

c at r or oz r 

1 dH \ d 

' -^irfrzJ) + ^if.zJ) = -{no + <yu-)H,+WzVo + ^zJ, (86) 



and 



c dt r dr dz 
1 dH^ Id d 

+ - ^{rfr<t,J) + -^{fz4,J) = -(ko + CTtr)i?0 + U'0?7O +^0-^, (87) 

c at r or oz 
where the individual components of S and ^ are given by eqs. (|33|l and (|39|) . 

4. MIXED-FRAME RADIATION EQUATIONS IN SPHERICAL GEOMETRY 

The spherical coordinates we use are the standard r, O, and 0. Here, <d is the polar angle measured from the z-axis 

and (j) the azimuthal angle. The unit vector in direction n is 

n — cos 9 Gr + sin 9 cos ?A ee 4- sin 9 sin -0 . (88) 

is the polar angle measured from the r-direction, and ij) is the local azimuthal angle between the projection of n 

onto the plane perpendicular to vector r measured counterclockwise from ee. The conservative form of the transfer 

equation in spherical coordinates is then: 

1 dl cos 9 d , ry , sin 6* cos d , . ^ ^. 1 d , . ^ „ . 

-^ + ^T- 0^ + ^-FTl^ sme/ — — (smHl) 

c at r^ or ^ ' r sm B aB r sm 6^ w ^ ' 

sin 9 cot Be) 

— (sm^/) = roo 4- roi 4- rio + rii , (89) 

r dtp 

where the individual right-hand side terms are given by eqs. H23|l - H26|l . specified for spherical coordinates, i.e. with 

rij — cos 9 Hr + sin 9 cos tp Hq + sin 9 sin tp H^ , (90) 
fijVj-' = cos 9v]r + sin 9 cos tp tVQ 4- sin 9 sin tpvu^ , (91) 

and 

UjtUkK^'' = COs9{lUrKrr + tVQKj-e 4- tJU^Krtj,) 

+ sin 9 costp{vjrKre + tveKee + w^if e<#.) 

4- sin 9 smtjj{tVrKr4, + tveK^e + lu^K^^) . (92) 
The corresponding moment equations read 

c Ot or r sm B oB 

= T]l^' - KqJ + ErHr + EqHq + E^H^ , (93) 

IdHr 1 d{r^frrJ) 1 



c dt r^ dr r sin B 9B 

frr 1 



{sin Qf re J) 



r 



J = -{na + <7ti )Hr + tUrVO + irJ , (94) 



IdHe ^ 1 djr'freJ) ^ 1 d 
c at or ?■ sm B oB 

/re-cotB/00 Mj , ~ I c T ^n^^ 

4 J = -(ko + tTtr)-ffe 4- we^ 4- Ce-/, (95) 

r 

and 

19i/^ , 1 dir^fr^J) , 1 9 



c dt r^ dr r sin B dQ 

fr4, 4- cot B /e0 



(sinB/^e^/) 



J = -(kq + o-tr)ff0 4- W0?7o 4- C^"^ • (96) 
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5. THE SPHERICALLY-SYMMETRIC CASE; EQUATIONS AND SOLUTION TECHNIQUES 

For initial tests of the scheme, we focus on the sphericaUy symmetric case. We introduce the usual notation /i = cos 9. 
We also assume that the velocity has only the non-zero component, Wr, which we denote as w. Due to symmetry, the 
only non-vanishing components of the vector H and tensors K and / are iJj., K„, and frr\ we denote them here as i?, 
K, and /, respectively. The transfer equation in the conservative form is written as 

ldl_ M djr'^I) 1 d[{l - ;u^)J] 
c dt dr r 



5^ — w [l + 5 — 



+0-0 

d\nH 
d\nv 



1 



dr r 

' [kq +<7o~ Hw(kq + (To 

9 In (To 9 In J 



2-5- 



SwiJ^ I 3 — 



91ni^ 
d In (To 



dlnH 



J 



dlnK 



d In i 



H 



K . 



The moment equations read 



IdJ 

c 9t 



1 d{r^H) 



and 



1 OH 
~c~dt 



1 d{r^fJ) 



dr 
1-/ 



no 



kqJ + :iH , 



(97) 



(98) 



(99) 



dr r 
where 5 = 2^ and i = ^r- 

We rewrite eq. (I97II in a non-conservative form, and express the time derivative through backward time differencing, 
while retaining the spatial derivatives, viz. 

1 - di 



dl 
or 



= -(Xo - fJ-Xi)! + iVa + mi ) + 



r dj-i 



where 



Xo = + cTo + 
Xi = ('^o 



1 



T 



r]( = (Tow 



no 



T 



-S 



cAt ' 

' cAi' 
3wrjo , 

= CTo , 

d In (To 



ainj 



-(Tow 



9 In d\nv 
dlriH 



1 + S- 



ao6w 3 



Vi = CTo<5 , 
9 In (To 



ainH 



and 



= (To(5u' 



9 In 9 in ly 
d\nK 



dim 



(100) 

(101) 
(102) 

(103) 

(104) 
(105) 

(106) 

(107) 
(108) 
(109) 

(110) 



Here, we split the appropriate coefficients into several parts depending upon which power of fj, it is associated with; 
subscript corresponds to the /^-independent part, subscript 1 to a linear term in fj,, etc. The superscripts J, H, K 
refer to the parts of emission coefficients that contain moments J, H, and K, and T refers to the "thermal" emission 
coefficient. 

In the tangent-ray approach, we consider transfer along the ray specified by a constant impact parameter, p. The 
coordinate along p is called s (we do not use the usual notation z to avoid confusion with z-coordinate used in 2-D 
cylindrical geometry), where 



,2n1/2 



(111) 
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Because the differential operator ^{d/dr) + r^^(l — fi^)d/dfj, is identically d/ds, one can integrate along straight lines 
because the characteristics of eq. (I1()0|I are straight lines. 

Because of the symmetry of the problem we consider only positive s. We denote the intensity propagating in the 
direction of increasing s by /"*" or /(/x), and that for decreasing s by /~ or The transfer equation along the 

tangent ray in the positive direction (increasing s) is given by 

QJ+ 

-g^ = -{xo - MxO-f"^ + {vl + ml) + 
ivo + mi) J + {v" + + f^'v2)H + miK , (112) 

while for the negative direction (decreasing s) it reads 

QJ- 

" -Q^ = -(Xo + fJ-Xi)!' + iVo - ml) + 

ivo - mi) J + ivo - mi + ^^\?)H - miK . (ii3) 

The transfer equation along the ray may be written in a compact form: 

^^^^ -IHp,^)- iP, s) , (114) 
X [P,s) ds 

where the corresponding total opacities and source functions 5*^ easily follow from the above expressions. The 
only complication is in the form of the source function, which may be written as 

S^{p,s)^S^ip,s) + a^{p,s)J{r)+b^ip,s)H{r)+c^ip,s)K{r), (115) 

where r — y^p^ + s^, and 

X'^iP, s) = Xo{r) - f^Xi{r) = Xo{r) ~ -Xi{r) , (116) 



X (P, s) ^ xa{r) + ^Xi{r) ^ Xo{r) + -Xi{r) , (117) 

snP,s) = ^ liM+lfZzM, (US) 

X^{p,s) X^{P,s) 

^-(p, ,) = !lL_J^ ^ vIir)- is/r)^T ^^^g^ 
X [P,s) X {P,s) 

+ ( ^ Vo + mi Viir) + is/r)vi , . 

"^"'^^^1^^^ X^P,s) ' ^'''^ 

and 

a- ip,s) = = ^i^^)-}'/l)^i , (121) 

X [P,s) X (P,s) 

and analogously for the other quantities entering eq. (|115|l . 

We introduce the following discretization. The radius grid is defined in terms of depth index d — 1, . . . , NR which 
increases inward from the surface: ri = i? > 7'2 > . . . > Rc, where Rc is the radius of the "core." The impact 
parameters are labeled by the same index as the radii, that is the impact parameter for the j-th ray is pj = rj. In 
addition, one introduces NC core rays with < Pnr+j ^ Raj — 1, • . • ,NC. The total number of rays (impact 
parameters) is, thus, NI = NR + NC. 

The moments J, i?, and K, which are integrals over angles, may be expressed as quadratures over the impact 
parameters, viz. 

NI 

J{rd) = J2 "^od [I^iPj, Sd) + I-{Pj,Sd)] , (122) 
j=d 

NI 

H{rd) ^^Wjdf^jd[I^{Pj,Sd) ~I^{pj,Sd)], (123) 

J=d 

and 

NI 

K{rd) = J2 ^^'^ [I^iPi^-^d) + I-ip„ sd)] . (124) 
j=d 

In the source function (eq. 1115(1 the parameters a^, b^, and are known functions of r, while the only 
unknowns are the moments J, H, and K which have to be solved for self-consistently with the transfer equation. Using 
the Feautrier scheme, one can in principle obtain an exact (non-iterative) solution, as in the case of the standard 
comoving-frame transfer equation in spherical geometry developed by Mihalas, Kunasz, & Hummer (1975). However, 
in the present mixed-frame approach with anisotropic scattering, the direct scheme is somewhat cumbersome. Even 
if one can solve the problem directly using the Feautrier formalism, it is nevertheless advantageous to use an iterative 
scheme. We shall, thus, outline an iteration scheme in 
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5.1. ALI iteration scheme 

The transfer eq. 1)114(1 . with the source function given by eq. (|115|) . is solved by an apphcation of an Accelerated 
Lambda Iteration (ALI) scheme. The solution of eq. H114() can be written formally as 

^ - E ^^^'St^' ' (125) 
d'=i 

where I^^ = I^{pj,Sd) and S'j^^ = S^{pj,Sd)- In other words, the specific intensity is understood as the result of an 
action of a certain operator (or matrix, when discretized) , A, on the source function. The iteration scheme adopted 
here is an application of the Jacobi preconditioning scheme, for which we write (dropping the superscripts ±): 

7-7 = A,- ,,5-7 + ^ HM'Sfi, . (126) 

d' = l 

In the usual astrophysical language, we use an approximate A operator given by the diagonal (local) part of the exact 
A. Using this expression, we can express the moments as 



Jr = ^ E ^.^Af^ddiS'^d + 4dJr + bfdHT"" + c%Kr-] + "old- terms , (127) 

NI 

= o E ^^d^'=d^tddiS% + ^tdJr + b%Hr + c%Kr] + "old" terms , (128) 



2 

3=d 



and 

r^ncw 



NI 

- ^ w,a^J:',dAtdd[S% + 4dJr + b%HT- + c%KT-] + "old" terms . (129) 



Here, the superscript ± indicates that we sum over both the + and — terms. After some algebra we obtain the set of 
three coupled equations for the new values of the three moments at the given radius r,: 



-A,« -Ki^< 
-Af l-A^ -A^^< 
-Ar -A^ l-A^ 



TIlCW 

•^d 






r JF 


- jr \ 






_ jjold ] 








(130) 


Tv'ncw 
^d 













where J^^ is given by eq. ((122|l with the specific intensity given by the "old" intensity 1°]^ (and analogously for H 
and if), and the matrix elements are given by 

NI 

Ai' = E ^^-^ HAtdd + ^jM^idd) > (131) 



NI 



Ai" = E ^^'i HdAtdd + bJ,dAj,dd) > (132) 

3=d 
NI 

Ai"" = E ^^'i (4Atdd + ^IdAldd) , (133) 

3=d 
NI 

A^' = E ^J.-i HdKdd - ^Aldd) ' (134) 

3=d 
NI 

^d" - E ^^'^ Kd^dd - bTd^^dd) - (135) 

3=d 
NI 

^d"" - E ^'^'^ /^^■.'i ^,d^.dd - '^IdA-dd) - (136) 



J=d 

and 



NI 

Af ' = E ^^■'^ ^'ld KiKdd + ^IdAj.dd) - (137) 

3=d 
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NI 

= E 4d iKAU + ^Id^ldd) ' (138) 

3=d 
NI 

^d"" = E ^^-^ 4d + ■ (139) 

To evaluate new values of the moments, one has to invert one simple 3x3 matrix per depth point. The individual 
matrix elements Aj^dd are evaluated during the formal solution step. 
The iterations procedure is as follows: 

(a) For given moments J^"-', H^'^\ .ftT'^"' (and with a suitable initial estimate of H^'^\ K^^^^), we perform a set of 
formal solutions for all impact parameters p, so we obtain new specific intensities, which we denote /pg(p, s). 

(b) We compute new values for the moments J^^, H^^,K^^ using eqs. p22|) - p24|) . with the specific intensity /pg. 

(c) We solve eq. H13()(l . radius by radius, to obtain new values for the three moments. 

(d) As long as the new moments differ from the old moments, we iterate steps (a) through (c) to convergence. 

Since the "acceleration" (step c) is very simple, the problem is effectively reduced to a set of formal solutions along 
the tangent rays. There are two main possibilities to solve the transfer equation along the ray, using either the Feautrier 
method, or the Discontinuous Finite Element (DFE) scheme. One can also use a 1-D short characteristics scheme; we 
have tested (i i7.1(l all three methods for the present study. 

We have found that one can use, without any deterioration of the ALI iteration procedure, a simplified preconditioner 
(approximate A matrix) that is obtained by dropping the off-diagonal terms of the 3x3 matrix A^^, namely by setting 

Ai" = Ar = Ai^ ^ Ar ^ Af ^' = Ar = O , (140) 

5.1.1. Tri-diagonal Operator 

A better preconditioner than that based on a local approximation of the exact transport operator is obtained by 
considering not only the local components of the operator A, but also terms corresponding to the nearest neighbors. 
In this case we write for the new moments 

1 

Jd ~ 2^ [ 'i'd.-l [^j.d-1 + "■j,d-vJd-l + OjA-l^d-l + S.d-l^d-l ) + 

3=d 

^td [Sla + + hf.HT'" + ctaKT"") + 

1 + 4,^,JT+i + + + "old terms" , (141) 



A± 



and analogously for and K^"^"". Generalizing the procedure described in ii5.ll we obtain a block-tridiagonal 

system (in the physical space) for the components of the moments; each block is a 3 x 3 matrix that couples all three 
moments. The diagonal block is the same as the approximate A matrix considered in ti5.ll and the off-diagonal blocks 
easily follow from eq. (|141|l . 

However, following our finding that the off-diagonal terms corresponding to the coupling of moments can be dropped, 
we end up with three separate tridiagonal systems in the physical space for the three moments J, H, and K. 

5.2. Augmentation of the ALI scheme by GMRES 

The ALI scheme outlined above can be significantly augmented by an application of a suitable Krylov subspace 
method, for instance the GMRES (Generalized Minimum Residual) method. There are a number of variants of the 
scheme; we have implemented the following one. Let us define the vector x composed of triads J, i?, K at all radii. Its 
dimension is, thus, 3 x NR. The general problem can be formulated as a linear system Ax — b, where matrix A is a 
block matrix of NR x NR blocks, each block being a 3 x 3 matrix, analogous to the matrix of eq. H130|l (containing 
also the off-diagonal elements of the A matrix). We define the "preconditioned residuum" at the n-th iteration, R^'^\ 
as a vector composed of — _ TjC")^ and if — if (") for all radii (i.e., a collection of the solution 

vectors of eq. (|130|) for all d). The GMRES scheme consists of consecutively finding "search vectors," P^'^\ whose 
products with matrix A are made orthogonal to the subspace spanned by the previously constructed search vectors, 
and which give the new estimate of the solution. 

The adopted algorithm is as follows: We start with x^"^ and set — i?'°^ Then, for each i = 0, 1, . . . , we compute 

(ApW,i?W) ^ ^ 
^ ^ (142) 



' (APW,APW)' 
x'''+^'^ = x^'^ + a,P^^ , (143) 
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R^'+'^^ = - a, AP^'^ , (144) 



h'ij — ^ 

and 



Pii=-—, A 7-rr-, J = 0,1,..., 2, (145) 



p(^+l) ^ ^ p^^pU) ^^4g) 

In practice, one can keep adding newer and newer search vectors. However, this may be cumbersome or too memory 
consuming. If so, one may actually stop and restart the orthogonalization process, or limit the orthogonalization to 
the k most recent search vectors. In this case, the summation in eq. (|145|l is replaced by X]}=i-fc+i PijP'^''^ ■ Such a 
method is sometimes called ORTHOMIN(k) (Klein et al., 1989). 

An important point is that one can (and should!) accomplish the procedure defined by eqs. H142() - (|146|l without 
performing explicit multiplications with matrix A, which in fact is never even assembled explicitly. It turns out that 
one can write 

= - , (147) 

and 

i 

J^p(^+l) ^ ^7^(^+1) + p^jAP^j^ , (148) 

so indeed explicit multiplications with matrix A are not needed. 

5.3. Discontinuous Finite Element (DFE) Scheme 

In 1-D, our implementation of the DFE scheme is a straightforward application of a method developed by Castor, 
Dykema, & Klein (1992), where the reader is referred for additional details and derivations. We present here the final 
formulae for the specific intensities, together with a description of how the elements of the approximate A operator 
are evaluated in the context of the DFE approach. 

In the direction of propagation we have recurrence relations for the finite elements: 

^d+i = o.d+i/2 (2 + Arrf+i/25d + bd+i/2Sd+i) , (149) 

and 

= O-d+l/2 + bd+l/2Sd - ATd+i/jS'd+l) , (150) 

where ^ 

ad+i/2 = + 2Ard+i/2 + 2) , (151) 

bd+i/2 = Ard+i/2(Ard+i/2 + 1) , (152) 

and 

Q+1/2 = 2(Ard+i/2 + 1) . (153) 
The specific intensity at point d is given as a linear combination of the discontinuous intensities, 

T Id^Td+1/2 + /^ATrf_i/2 

It can be shown (Castor et al. 1992) that this choice of linear combination of the discontinuous intensities makes the 
scheme second-order accurate. 

As follows from the above expressions, the diagonal and first off-diagonal elements of the transport matrix A are 
given by linear combinations as in eq. H154|) : 

Ad,,Arrf+i/2+Atj-ATrf_i/2 . 1 ,,,,, 

Ad,j = — ^^-r — r-^^^ , J = d-l,d,d+l, (155) 

where 

K+i,d+i = ad+i/2&d+i/2 , (156) 



K+i,d = «rf+i/2 (ATrf+1/2 + 2A^_^ j , (157) 

^d,d ^ ad+i/2ibd+i/2 + Cd+i/2Kd) ' (158) 

and 

^d,d+i = -A^d+i/2ad+i/2 • (159) 
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5.4. Feautrier scheme 

While using the short-characteristics or the DFE schemes represents a straightforward appHcation of these methods 
(for a review, see, e.g., Hubeny 2003), an apphcation of the Feautrier scheme requires a generahzation of the Mihalas, 
Kunasz, & Hummer (1975) formahsm, which we describe in this section. 

To use the Feautrier scheme to solve the transfer equation, we introduce the usual symmetric (mean-intensity-like) 
and antisymmetric (flux-like) Feautrier variables 

C/ = i (/+ + /-), (160) 

and 

1/ =!(/+-/-). (161) 

Eqs. (|112|l and (|113|l can then be written as a set of two differential equations for U and V: 

dV 

— = -xoU + i^XiV + r,l+ + r,oO + (r,o^ + H , (162) 

and 

— = -x^V + nxiU + Til- + M ('7? + riiJ + ViH + if) , (163) 

where 

%^+ = %^ + ^[/°(/^)+/°(-/i)] , (164) 

and 

v^- ■ (165) 

Here, we take into account the fact that tjq depends on /i only through the intensity at the previous time step, 
The symmetric and antisymmetric averages are, thus, denoted as rj^ ^ and rj^ ~, respectively. 
Furthermore, we introduce a modified optical depth along the tangent ray, 

dT = —Xods , (166) 

and the "source terms" S'^ and S~ , where 

S+ = S^+ + S^J+{So + ^i^S2)H, (167) 
5- = S^" + niSf + 5f J + 5f iJ + S^K) , (168) 

and 

5*^ = ^, for i = 0,l,2 and X=^T,J,H,K. (169) 
Xo 

The equations for U and V can then be written simply as 

dV 

— = U -aV - S+ , (170) 

OT 

and 

'~\TT 

— = V -aU - S- , (171) 

OT 

where 

a = Ai — . (172) 

Xo 

To reflect the mean-energy character of U and the flux-like character of V, we stagger the U and V meshes by half 
a zone. That is, U is taken with integer indices d, while V is taken with half-integer indices d ± 1/2. The discretized 
eqs. (|T7n)l and ((T7T|l read 

- Ud - -z-{yd+i/2 + Vd-1/2) ~ i>d ' (1''^) 



and 



where 



ATd 2 

- - ^iUd,^ + Ud) - S-^,,, , (174) 

= \{sd - Sd+i)[{xo)d+i + (xo)d] , (175) 
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(176) 

(177) 



The half- integer /irf+1/2 is given by an equation analogous to eq. (|177|l . where the corresponding radial points are given 

by 

1 



(178) 



That is, the cell center is defined in such a way that the volume at is half of the volume between Vd and r^^i. 

Discretized equations may be written in a more compact way by introducing 



1 



ATd 2 ATd 2(x. 



1 ^ Md(xi)d 



d 



(179) 



where d can assume integer or half-integer values. Equation (|174|l can then be solved for Vd+1/2 to read 

Vd+1/2 = Ud+ib+^^^^ - Udh^^^f^ + 5*^+1/2 • (180) 

Analogously for V"d_i/2, we have 

(181) 



d-1/2 



Using these equations, we can eliminate Vd±ii2 from eq. (|173|l to obtain 



Ud-ihd_i/2bd + Ud[l + bj_^/2bd 



"d+l/2"d 



^d + '^<i+l/2 



^d ^d-1/2 ■ 



(182) 

The inner and outer boundary conditions can be handled in straightforward fashion. For the inner boundary 
condition, d — D, we have /+ — I~ and /i = 0. We, thus, have V = and S~ = 0, the latter equality following from 
eqs. ()168|l and (|165|1 . The two Feautrier equations are written 



D 



dT 



f 1 



D 



and, thus, 



9r2 



Ud-S} 



dS_ 



dT 



The 2nd-order form of the boundary condition follows from a Taylor expansion of Ud-i around D, 



dU 

Ud-i = Ud - Atu-i/2 [ 



2^'D-1/2\q^2 



Using eqs. (|183f) - (|185ll . and expressing dS /Ot as a difference, we obtain 



U. 



D 



1 



'_D-l/2 ^ 



Ud-1 



S 



1 



D 



' D-1/2 



Atd- 



S 



1/2, 



D-1 ■ 



(183) 
(184) 

(185) 
(186) 

(187) 



'D-1/2 

For the outer boundary condition at d = 1, we take = 0, and, thus, Ui — Vi. In order to write the 2nd-order 
form of the boundary condition, we first write a general 2nd-order equation that is derived from eqs. (|170|) and H171|l : 

(188) 



9r2 



iS — 



dT 



The 2nd-order form of the boundary condition follows from the Taylor expansion of U around d = 1, 



U2^Ui + AT3/2 — 

Substituting eqs. (|170|l and (|188|l into eq. (|189|l . we obtain 



1 



A-2 



92 [/ 



da 2(1 -ai) 



Ar3/2 



At2 
^^3/2 _ 



3/2 Q^2 



-U2 



(189) 



^^3/2 



^'^3/2 



ai ) 



dS_ 



(190) 
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Finally, expanding the derivatives of a and S we obtain 



AT3/2 AT3/2 A 



3/2 



U: 



2 



^^3/2 



St + ^r^ ~a^S^. (191) 

Equations (|182|l . (|187|l . and p91|) form a tridiagonal system that is solved by the standard elimination method. 

An alternative, and in fact more accurate, way of formulating the Feautrier scheme is by introducing the integration 
factor, q, defined by 

-^ = a, (192) 
q OT 

with which the Feautrier eqs. H17U|I and H171|) are rewritten as 



1 d{qV) 
q dr 

and 

1 d{qU) 



= U-S+, (193) 



^V-S~ . (194) 



q dr 

Using these equations, we can eliminate V, to end up with a single second-order equation for U: 

=lU- qS+ ^ . (195) 

This equation is discretized in the standard way, namely for d = 2, . . . , Z? — 1 we have 

_ qd-iUd-i qdUd f 1 1 \ _ qd+iUd+i ^ ^+ gd+i'S'^+i - qd-iS^_^ 

Ar<jAT<j_i/2 Ard ^Ar<j_i/2 At^+.^J ATdATa+1/2 '^'^ ^Au ' ^ ' 

where we represented the derivative d{qS^) / dr by a centered difference. 
The discretized boundary conditions are obtained similarly as above, and are given by 



1^Ui\l + -r— + ^rj2-\-q2U2^rj^ = qiS^ + -r , (197) 

Ar3/2 At^^ AT3/2 A3/2 



and 



qnUo ( 1 + ) - <1d-iUd-i-^^ - qoSi 2£-lEl^ . (198) 

5.5. Moment equations 

At first sight it may seem unnecessary to deal with the moment equations because all the necessary information 
about the radiation field is provied by the specific intensity, whose evaluation was described in the previous parts of 
^ However, there are several reasons for considering the moment equations: 

i) The hydrodynamical equations, and in particular the energy and momentum balance equations, and the electron 
fraction equation (for core-collapse supernovae), contain only moments of the specific intensity, not the intensity itself. 
It is, thus, natural to work in terms of the radiation moments. 

ii) It is actually more consistent to work in terms of moments, since in the moment equations the angular integrations 
are performed analytically, while when using the specific intensities the radiation moments are evaluated by a numerical 
integration. Generally, these two ways give somewhat different results, and one should carefully assure the consistency 
between the radiation moments used in the formal solution of the transfer equation, and those used in the material 
equations. 

iii) There is a practical aspect. Solving moment equations is considerably faster than solving the full angle-dependent 
transfer equation. In fact, the full angle-dependent transfer solver may be viewed as a tool to provide just the Eddington 
factor, while the moments are obtained by a numerical solution of the moment equation. In many cases the Eddington 
factor changes only slowly with time, and therefore one does not actually have to update the Eddington factor in 
each timestep; instead the Eddington factor may be held fixed for several timesteps. This may obviously lead to a 
significant reduction of the computer time required. We shall return to this point in ii7.2l 

The moment equations in conservative form are written as 

1 dJ 1 d(r^H) , , 

c ot or 



and 



where 



and 
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c dt dr 



KO + (Jtr 



■J = -(ko + crtr)i^ + Wl-jo + , 



9 In (To 5 In 
9 In 9 In J/ 



— (2-5- ^ ^° - + /kq + /(To (1 + ^ ^'^ ^° + 5 d\nK 

3 V 9 In J/ dhiv 1 \ dlnv ?> d\nv 



w . 



We rewrite these equations in a more compact and useful form, again expressing the time derivative throu; 
time differencing: 

dh 



— = kj - (3h- Sj , 



and 

where the modified moments are given by 
and 

and where the auxihary quantities are given by 



dr 



= h - aj - Sh 
h = r^i? , 



gJ ^ ILL s" ^ — 



KH 



KH 

2 th I ' -"0 2 ~ , ''^-ffo 

/ 1 \ 1 1 



^ e + (!-/) A ^ s 

a = , P= , k 

JHH KH KH 



cAt ' 



and 



Discretization of eqs. H203|) and 1204() yields 

hd+1/2 — hd-1/2 



dr = —Kndr . 



and 



jd+i - jd 



= kdjd - ^{hd+1/2 + hd-1/2) - •S'i . 



h-d+i 



/2 



2 

"d+1/2 



O'd+l + jd) - "5*^+1/2 ■ 



ATd+1/2 

We eliminate h from eq. (|213|l to obtain 

/»d+l/2 = Jd+lld+1/2 - 3dla+l/2 + '^H 

Performing the same operation for h^_i/2 and substituting into eq. (|212|l . we obtain 

- Jd-i Cdld-1/2 + jd [kd + 47d+i/2 + 7^-1/2) - jd+i chd+1/2 



^ Cd^H ^ '^d '^H ' 



where 



and 



1 0d 
ATd 2 ' 



± 1 , ad 
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(200) 
(201) 

(202) 
;h backward 

(203) 

(204) 

(205) 
(206) 

(207) 
(208) 
(209) 
(210) 
(211) 

(212) 
(213) 

(214) 



(215) 
(216) 

(217) 



and where we understand that the depth indices in eq. H217|l have half-integer values. 
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The boundary conditions are expressed using the second-order form of the moment equation. EUminating h from 
eq. (|2()4|l and using eq. (|2()3|l . we obtain 



k ~ ^ + a'^ ] j ~ {a + f3)h 



Sj - aS_ 



dS 



H 



H 



dr 



The outer boundary condition (at d = I) follows from expanding j in a Taylor series: 



•1 A ^-^ 
j2 = Ji + At — 



1 A 2 d^j 
-At^ — - 

2 dr2 



(218) 



(219) 



Substituting for derivatives from eqs. H203|l and H218() . we obtain 



Ji 



A: 



Ih 



( 2 



3/2 



■A 



AT-ii2 



ki + al 



' J2 



(220) 



where we have also used the 1-st order form of the derivatives of a and Sri [for instance, {da/dT)i = {a2 — Q!l)/AT3/2]■ 
Finally, 

fH^^^^, (221) 
Ji JiJi 

is the "flux Eddington factor" which is evaluated in the full (angle-dependent) transfer solution. 

The inner boundary condition (d — D) is analogous. There, we use eq. (|218f) and the symmetry condition hn — 
to obtain 



JD 



^^D-l/2 -^^15-1/2 



]D~1- 



' D-1/2 



qD I qD-1 

bj - ao^H - 



(222) 



Atd_i/2 

Equation (|215|l . together with boundary eqs. H22U|1 and (|222|) . form a tridiagonal system which is solved by a standard 
Gaussian elimination (also called a forward-backward sweep). 

5.5.1. Matrix representation 

To prepare for the formalism used in !|7| we note that we can write down the tridiagonal system for j defined by 
eqs. (|215|l . (|220|l . and H222(l in a matrix form 

T ■i = S' + U -S" , (223) 

where j = (ji, J2, ■ ■ ■ , ju)"^ is a column vector of moments jd at all depth points. There is an analogous expression for 
the source function vectors S"^ and . Matrices T and U are tridiagonal matrices. The matrix elements of T are 
given by 

= -c^7<ili/2 - ^'d,<i = ^<i + 47d+i/2 + CdTjLi/2, TdM+1 = -cj-fj+y^ ^ (224) 

for 2 < d < I? — 1. Expressions for d ~ 1 and d ~ D easily follow from eqs. (|220() and H222() . The matrix elements of 
U are given by 



Ud.d-i 



Un = 



2 
1 



AT3/2 



2 ' 

' ai 



Ud.d+i 



U12 = 



2 < d< D -1, 



AT3/2 



(225) 
(226) 



and 



Udd 



ao , Ud.d-1 



D/2 



(227) 
(228) 



Ar£i_i/2 

One can formally write a solution for j as 

j = ■ S-' + T-^U -S" = Ao- S'^ + Al • , 

with new matrices Aq = and Ai = T~^U. 
We can write a similar equation for vector h — {hi,h2, ■ ■ ■ ,/iD)^ using eqs. H214|) . H228|l . and expressions hd = 

{hd-1/2 + hd+1/2)/'^, hi = /hJi, and hd = 0: 

h^V-i + WS" ^ {VAo) ■ + {VAi + W)-S" , (229) 
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where V and W are tridiagonal matrices whose elements are given by 

F.,.-r = -%^, ^.,.+1 = %^, 2<d<D-l, (230) 

Vii = fn , Vi2 = Vdd = Vd,d-i = . (231) 

Wdd = l Wd,d±i = \, d<D-l, (232) 

= Wdd = . (233) 



and 



We write the expression for h as 



h = A2 • S'^ + A3 • , (234) 



where A2 = VAq and A3 = V^Ai + 14^. 
Finally, the true moments J and H are given as 

J = diag(l/r2/rf).j, i.e. Jd = ^ , (235) 

and 

H = diag(l/r2).h, i.e. Ha^r^^ha. (236) 

5.5.2. Generalized Sphericity Factors 

It turns out that the discretization of moment equations represented by eqs. (|215f) . H22U|I and (|222|) . may lead to 
significant inaccuracies in the evaluated moments j and h. This is essentially due to a first-order representation of the 
Ph and aj terms. An improved numerical technique easily follows from introducing integration factors of eqs. (|203|l 
and (|204() . In the classical photon transport, such an approach was pioneered by Auer (1971), who coined the term 
"sphericity factors." 



We rewrite eqs. \M6\ and (|204|l as 



1 d{qhh) 

1 ^ kj - Sj , (237) 

qh dr 



and 



where the integration factors are defined by 



and 



1 d{qjj 



qj dr 

djln qh) 
dr 

d{lnqj) 



- Sh , (238) 
P (239) 



dr 

We use eq. (|237|l to eliminate h from differenced eq. (|238|l . which leads to a second-order equation for qjj: 



dx^ qj 



— H 1 

qj dx 



(240) 



(241) 



where 

dx = ^dT = -^KHdr, (242) 
Qh Qh 

represents a modified optical depth increment. Equation H241() is second-order accurate. Discretization of this equation, 
and the second-order boundary conditions, are quite analogous to the classical Feautrier scheme. 

6. IMPLICIT COUPLING OF RADIATION TO MATTER 

We outline here a procedure for the implicit coupling of radiation and matter that is based on an application of the 
Accelerated Lambda Iteration (ALI). The procedure is a generalization of a treatment developed by Burrows et al. 
(2000); an analogous procedure was developed in the context of static stellar atmospheres by Hubeny & Lanz (1995), 
who coined it the "hybrid CL/ALI" method (CL stands for Complete Linearization). 
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6.1. ALI treatment of radiation 

The specific intensity can be written as 

= Ap.^p. , (243) 

which is just another way of expressing eq. H125|l . Here, we do not discuss energy redistribution and inelastic scattering, 
but provide our formahsm for the coupling of energy groups in the Appendix. With independent energy groups, we 
drop the energy subscript, u, but remember that all quantities still depend on energy. 

We assume that at the end of a given time step we know the source function and the corresponding specific intensity, 
which we denote as "old" quantities, I"'*^ = A^5'°''*. We express the "new" specific intensity as 

/-- = A^(5°i<i + ^5,,)- (244) 

We wiU also denote I^''" = /"'"^ + (5/^, and express A^ = A* + 5A^. Substituting these into eq. H244|l . and neglecting 
the 2-nd order term SS SI we obtain 

SI^ = a; SS^ . (245) 

This states that the correction to the specific intensity is given by the action of an approximate operator on the 

correction to the source function. 

The source function is given by eqs. I|167|) and Hl()8|) . which we rewrite here as 

S^^So + a^J + b^H + c^K , (246) 

where 5*0,0^,6^, and are known functions of the material state parameters T, Y^, and p. The correction to the 
source function can, thus, be written 



where, for instance. 



'^•^M = + -^SYe + -^5p + a^5J + h^6H + c^5K , (247) 



9T - 9r + 9r + 9r ^ + ar ^ ' ^^^^^ 

and analogously for the other derivatives. The expression for 51^ easily follows from eqs. H245|l and H248(l . Integrating 
over angles, we obtain equations for the corrections to the moments (5 J, 5H , and 5K . For instance, for 5 J we obtain 

+5J j A*^af,dp/2 + SH J K*^b^dii/2 + 5K J A;c^d^/2, (249) 

and analogously for 5H and 6K, the only difference being that the angular integrals are modified to / ■ • • pdfj, and 
/ • • • /i^d/i, respectively. 

The simplest choice for the A* operator is the diagonal (local) part of the exact A, as discussed in ^5.11 Notice 
that Burrows et al. (2000) used a tridiagonal part of the exact A, which generally gives faster convergence. A great 
advantage of the diagonal operator is that it yields relatively simple expressions in 1-D, and is very easily generalized 
to multi-D, where anything but a diagonal operator would lead to somewhat cumbersome expressions that involve a 
number of neighboring cells. We will consider both cases, a diagonal as well as a tridiagonal approximate operator. 

With the choice of diagonal operator, eq. H249|) . and analogous equations for SH and 6K, can be written as separate 
systems for each spatial zone (depth point): 

/5Jd\ 

\skJ 

■ST+\ E,^-^ \ SY,+ \ Y."^" \-5p , (250) 

where we observe that the elements of the moment-coupling matrix are exactly those defined by eqs. H131|) - (|139|l 
needed in the formal solution of the transfer equation, and where the elements of the vectors representing the appro- 
priately angle-averaged derivatives of the source function are given by 





K =2] /m^^'^/^' (251) 



^1-" = ^/'^;^^^, (252) 
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2 A * ^"^M J 



(253) 



and analogously for T,^"'"^ , etc. The essential point here is that we are able to express 5 J {SH, 5K) as functions of 5T, 
SYf,, and Sp, through a simple 3x3 matrix inversion. 

These equations are quite general. In practice, we perform an implicit update for the temperature and Ye only, 
because the density is most naturally updated in the operator-split fashion in the explicit hydro step. In this case, we 
formally set Sp = 0, and write (dropping the depth index d): 

(and analogously for 6H and SK) where 

7" 1 

There are analogous expressions for 5'^'"', etc. 

There is an alternative way of formulating the ALI-based corrections dJ and SH, namely those based directly on 
the moment equations H228|l and (|229|l or H234() . together with eqs. (|235|) and H236() . Since the material equations are 
written in a form completely consistent with these moment equations, this approach seems to be more suitable. We 
write for the correction SJd, again using for the approximate operator the diagonal (local) part of the exact operator: 



-Ar 




^d 



(254) 



(255) 



: fd SJd 



(Ao; 



dd 



OTd + -ITTT- OYed 



dT " dY, 
Thus, the auxiliary quantities 5* can be written 

,T,J 1 



(A 



{^o)dd5Si 



H 



l)dd 



dT 



STd 



{Ai)dd SSd 



dY, 



(Ao) 



dd 



and 



^dfd 



iM)dd 



dS-' 

f ds-' 

\dY, 



(Ai) 



dd 



(Ai) 



f ds" 

\ dT 



dd 



dY, 



(256) 



(257) 



(258) 



Here, the integrations over angle have been performed analytically, as opposed to the numerical integration of the 
previous approach. 



Analogously, for 5H one obtains 

fdS-' 



(A2) 



dd 



V dT 



STd 



dY, 



rlSHd = {P^2)dd5Si 

H 



SYe. 



(1^3) dd 



dS 



dT 



STd 



{A3)ddSS!^ = 
dS 



H 



dY, 



■SYe. 



and, thus, 



and 



T.H 



(A 



2}dd 



{M)dd 



fdS-^ 



\ dT 

f dS-' 
\dY, 



(A 



3)dd 



dS 



H 



dT 



(A3) 



dd 



\dY, 



(259) 



(260) 



(261) 



In order to evaluate these coefficients, one needs first to compute the diagonal and a few off-diagonal elements of the 
inverse of the tridiagonal matrix T. This is done very efficiently using the procedure suggested by Rybicki & Hummer 
(1991). These elements are evaluated during the solution of the original set of moment equations and comes at almost 



no additional cost. We, thus, take elements T^^d , T^ and T^ ^_|_2 as known, and write 



{Ao)dd = 



(Al)d<J : 

(A2)dd 



Td,d-iUd-i,d 



dd 

T,,,, Udd 



'-dd 

VddT 



dd 



Td.d+iUd+i.d 



Vd^d+iT 



-1 



d+l.d 



and 



{^3)dd — Vdd (,^l)dd 

Vd^d+1 



Vd,d-i {Tdli,d-iUd-i..d + Tdli,dUddj + 



I rp-l 

\ d+l. 



dUdd 



T 



d+l,d+l^d+l,d 



dd ■ 



(262) 
(263) 

(264) 



(265) 



Numerical experience has shown that while the first way of evaluating the moment "derivatives" ^ is fast and simple, 
the second way is more consistent and accurate. This is because the material equations (energy balance and electron 
fraction equations) are based on the equations for the radiation moments, in which the integrations over angles are 
done analytically. Thus, the approximate operator based on the radiation moment equations can naturally be used 
in the material equations, because in the first way of computing the approximate operator the necessary angular 
integrations are performed numerically. 
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6.2. Linearization of the energy and the electron fraction equations 
As follows from the analysis presented in ^2A\ (eas. 1491 and l^n|l . the energy balance equation is written as 
DT f°° 

pCv— = -inJ2j^ [nl{v,T,Y,)-Kl,{v,T,Y,)Jl + Xr{v,T,Y,)Hl\dv, (266) 



where the summation extends over all neutrino species i {tjq is the corresponding rj^^), and where 



Xt = w [ 2ko + ° ) = w (ko + Ko) , (267) 



The electron fraction equation is 
DY 



rNa 



Dt 



= 4^ V / [i^l^iv, T, Y,) ~ ^Uv, T, Y,)Jl + X^{v, T, Y,)Hl] — , (268) 
"'0 ^ 



where — —1 for z^g neutrinos, = 1 for i^e neutrinos, and s' = for other neutrino species, Nj^ is the Avogadro's 
number, and 

Xy^wko- (269) 
The implicit (backward time differencing) forms of these equations are 

r Hii^,T,Y.^~^hi^,T,Ye)Jl + XU'^,T,Ye)Hl]diy + Am{T)^ 0, (270) 

and 

^Z'^^^aI" - E [^o(^, T, Y,) - nl^iv, T, Y,)Jl + X^{v, T, Y,)Hl\ ^ + ADV(ye) = . (271) 

Here, ADV(T) and ADV(ye) are the formal advection terms. Depending on the overall hydro scheme to which the 
present formalism is being implemented, the actual advection terms may be already considered in the hydro step, 
and not in the implicit update of T and Y^. due to radiation, which we consider here. In this case, we have to 
set ADV(T) — ADV(ye) = 0. For certain testing purposes, we may consider the advection terms as a part of the 
present implicit update. Therefore, we write the advection term at radial zone d (at radius r^), using a second-order 
representation of the spatial derivative, as: 

ADV(T)<j = Z [ed.d-iT{rd-i) - {ed.d-i + ed.d+i)T{rd) + ed,d+iT{rd+i)] , (272) 

47r 

and analogously for ADV(ye). Here, Z = if the advection was already treated in the hydro step, and Z = \ li 
the advection is being treated as a part of the implicit update. The coefficients e depend on an adopted form of 
differentiation formula. For a simple centered difference, 

ed,d+i = ~e.d.d~i = l/(''<i+i - rd-i) ■ (273) 
One can also use a truly 2nd-order formula, where 

ed.d+i - 7 77 r , ed.d-i = -7 77 7 • (.^74) 

{rd+i-rd)(rd+i-rd-i) [rd - rd-i)(rd+i - rd-ij 

They can be treated either explicitly or implicitly. If they are treated explicitly, and if we use a diagonal A* 
approximate operator to treat radiation quantities, then the problem can be formulated as a local problem in space. 
We shall first consider the local formulation, while we will consider the non-local formulation, with an implicit treatment 
of the advection term, later. 

We introduce the state vector 1/; = (T, Ye)^, and write the system of material equations H270f) and (|271|) as 

P(V')= 0, (275) 
which we linearize using the standard Newton-Raphson technique, 

P'(V'^"^) • Si:^"^ = -P(?/>(")) , (276) 

where P' is the Jacobian of the system; 

i.e., the i,j element of the Jacobian is the partial derivative of the i-th equation with respect to the j-th variable. 
In expressing the Jacobian, we use the expression SJ — {dJ/dT)ST + [dJ / dY^)6Yi,, together with eq. (|254|) . namely 
5J = ip'^'^^ST + ip^'-'^SYe, which allows us to associate dJ/dT ~ xp'^'-^ and dJ /dY^ — ^i^<"-^ . We can write analogous 
expressions for dH/dT and dH/dYe. The energy and electron fraction equations do not couple the adjacent depths. 
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and, thus, the Jacobian has a simple block-diagonal structure; each block is a 2 x 2 matrix which we denote A (we 
drop the index indicating the radial zone). The matrix elements are then given by 

^ ^ + E - f ^0 + ^i^o ~ ^0^-- + A,^-^') , (278) 

^21 = - E [ - #-^0 + ^i^o - -0*"^^''^ + AV^-'--) ^ , (280) 

and 

^^^~4^"^'X 91;^°"''°* jV' ^^^^^ 

The quantities that depend on T and Ye, such as k*, 77*, and S% are to be evaluated at the current iterate T^^"^ and 
yi"-* . The starting estimate is obviously given by r'"' = To and Y^°^ = 1;°, that is by the values at the end of previous 
time step. 

The individual blocks of the right-hand-side vector, P{ip^'^''), which we denote as &, are given by 

^ At ° " ^ / ^"^^ " ""-^^ + " ADV(To) , (282) 

i ^ 

and 

Notice that we have to treat here the advection terms (if they are taken into account, Z = 1) explicitly. 

The problem is, thus, reduced to solving, for each radial zone, the system ASij: = b. However, the moments, as well 
as their "derivatives," V?, are evaluated at the end of the previous time step. 

In fact, one can avoid inverting a 2 x 2 matrix A by setting its off-diagonal elements to zero: A12 = A21 = 0. 
We have verified that this does not decrease the convergence speed in any appreciable way. We can also introduce a 
slightly modified notation and write An = A22 — -B|^, bi = and 62 = Rd (where all the state parameters 
and radiation moments are taken at depth d), and write the global linearization scheme for the corrections as 

B^ST = , B^5Y^ = , (284) 

where 

{B'^)d,d' = Bj^5d,d' , {B^)d,d> — Bj^jSd^d' (285) 

are diagonal matrices and R^ = {Rf > R2 ' ■ ■ ■ > -^dI' ^^'^ — {^^i, ST2, . . . , 6Td} (and analogously for Ye) are the 
appropriate column vectors. 
The iterations proceed as follows: 

(a) Taking current estimates of J* and (and taking Jq and Hq as starting values), we perform the inner Newton- 
Raphson iteration loop, eq. I|276(l to determine estimates of T and Ye at the end of the given time step. 

(b) We perform a formal solution of the transfer equation using newly computed values of T and Y^ to compute new 
estimates of J' and H\ and possibly also new vE'^'J'. 

(c) Return to step (a) and iterate. This is an outer iteration loop. 

In this nested iteration loop, one can in principle obtain a fully self-consistent implicit solution for J, H, T, and 
Ye, i.e. for the radiation moments and the material quantities. However, this may be too time-consuming (although 
possible in 1-D), so we usually resort to approximate schemes. The most useful approximations (roughly in the order 
of decreasing overall importance, i.e. their expected influence in saving computer time) are: 

(1) In the outer iteration loop, one does not perform the full transfer solution for specific intensities. Instead one 
only updates the radiation moments, keeping the Eddington factor fixed. This means that the Eddington factor is 
going to be treated explicitly, while the rest of the quantities (radiation moments and material quantities) are treated 
implicitly. In this case, the Eddington factor may be updated at the end of the nested iteration loop (at the end of 
the given timestep). It may even be kept fixed for several time steps (which will almost certainly be necessary in the 
2-D case). 

(2) We may reduce the number of iterations in the outer loop, perhaps even to 1 (only recalculating the radiation 
moments after the first implicit update of T and Y^ is done. 

(3) In the inner loop, we may hold some quantities fixed. For instance, we may update rj, k, S, and their derivatives, 
but keep the 4* factors, that depend on the approximate A operator, fixed. 
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(4) The extreme variant of the above strategy is to use the so-called Kantorovich variant of the Newton-Raphson 
procedure that consists in keeping the Jacobian fixed altogether. Since setting up the Jacobian is the most time- 
consuming part of the inner loop (its inversion is easy, since it is only a 2 x 2 matrix, or just two divisions, if we use 
diagonalized matrices), this may lead to considerable savings. 

The essential point behind the two last simplifications is that one needs solve only the material equations H27()|l and 
H271|) (together with the transfer equation) exactly. Since the Jacobian is only one possible means to obtain the exact 
solution, there is no need to compute the Jacobian "exactly" from its mathematical definition. In other words, the 
basic linearization equation l|276l) computes only corrections to quantities, not quantities themselves, and, thus, one 
can afford approximations. Here, any matrix will do, as long as the process converges sufficiently quickly. 

6.3. Linearization using tridiagonal operator 

A useful generalization of the previous formalism is to replace the diagonal representation of A* by a tridiagonal 
one. Although it would be easy to retain the coupling of ST and 5Ye, which would lead to a block-tridiagonal system 
for the corrections, we consider here the case of uncoupled ST and SY^, as in eq. (|284l) . Moreover, we can now treat 
the advection terms implicitly at no cost, because the matrices already consider the radial zone coupling. Specifically, 
eq. H284(l remains valid, but matrices and become tridiagonal, where the diagonal elements are the same as 
before, and the off-diagonal elements are given by 



and 



where 



and 
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(287) 



(288) 



(289) 



and where X stands for T or Ye. Obviously, if we adopt the explicit treatment of the advection terms, the last terms 
of eqs. H286(l and H287(l disappear, and appear in the right- hand-side vectors i?'^ and . 

6.4. Treatment of the optically thick region 

Since the radial optical depths close to the center may be very large for most neutrino species, we must assure that 
we recover the diffusion limit exactly. This must apply not only for the solution of the moments themselves (which 
we are doing by construction since we are using methods which are second-order accurate), but also in the implicit 
update. To this end, we found it best to replace the above outlined approximate operator by that based on the diffusion 
approximation. Such an approach is often used in neutron transport theory where it is known by the name Diffusion 
Synthetic Acceleration (DSA). In other words, we employ the approximate operator, not as a diagonal or tridiagonal 
part of the exact transport operator, but rather as the corresponding expression in the diffusion approximation. 

In this case we use the approximation 



J 



where S = r^g /kq, and we set r 



Jd — Sd 



1 
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1 



1 d^S , 1 dS 

+ 3^' and H^--, 

In the discretized form (for 1 < d < NR.) 
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and (using a centered numerical differentiation) 
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(290) 



(291) 



(292) 



We stress that we do not use these expression for evaluating the moments themselves (this is done by exact solution 
of the moment equation) , but only to evaluate the derivatives dJ/ dT, dJ/ dVe (and analogous ones for H) to be used 
in the implicit update. In our notation, we write: 



^Vd = 



1 - 



1 



3 ATd \ATd 



1 



-1/2 



1 



1/2 



dSd 
dTd 



(293) 
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and 



^T,J I 9Sd±i 



and 

,T.H |_ 2 dSd±i 



<f±i = ±^^. <'f = 0. (295) 



and analogously for and "if^'^ . 

To assure a smooth transition between the previous formalism, which should be used in the optically thin regime, 
and the present one, we set the "derivatives" ^E" equal to a linear combination of the original quantities which we 
denote by ^'o, and the present ones, which we denote by ^'dsa (for Diffusion Synthetic Acceleration): 

= eM~r)^o'' + [1 - expi-r)]^l'^^ . (296) 

This choice of the linear combination is not unique, and some other choice might be better, but we found it to be quite 
robust. 

7. NUMERICAL TEST AND RESULTS 
7.1. Stationary solutions 

We first study the convergence pattern and the behavior of the solution in the stationary case. We consider a 
fictitious structure, based on an angle-averaged snapshot of a VULCAN/2D simulation taken from Dessart et al. 
(2006) at 200 ms after bounce. The original 2-D structure was appropriately averaged to yield a 1-D, spherically 
symmetric structure. Figure Q displays the temperature (T), density (p), and the electron fraction (Ye) as a function 
of radius. 

In the present tests, we consider three neutrino species, v^, v^, and a composite of Vr, v^, and 9^, which we 
denote as 'V^^." For each species, we consider 16 energy groups logarithmically equidistant between 1 and E'max, where 
Eraavi = 300 MeV for Vf. neutrinos, and iSmax = 100 MeV for the other two species. 

First, we test the global ALI iteration scheme. We consider five cases: diagonal and tridiagonal operator, with or 
without the GMRES augmentation, and another acceleration scheme related to the GMRES, namely Ng acceleration 
(Ng 1984; Auer 1991; Hubeny & Lanz 1992), which is widely used in astrophysical photon transport work. The solver 
for the formal solution of the transfer equation used in these tests is the DFE scheme. We stress that we consider here 
the most stringent test case of the ALI iteration scheme, the one for which the initial value of the radiation intensity is 
set to zero, and we compute the stationary solution directly. In a realistic time-dependent solution, one starts with the 
radiation intensities at the previous timestep, which are already relatively close to the solution at the given timestep, 
so this generally requires far fewer iterations. 

Figure m displays the convergence pattern for Vf. and neutrinos with E = 8.6 MeV, for the 5 different setups of 
the transport solver. Convergence is pretty fast for v^, neutrinos; an application of a tri-diagonal approximate operator 
leads to a somewhat faster convergence, as does the application of the GMRES scheme. The case of neutrinos is 
more interesting. Here, we see that the tridiagonal operator yields a substantial speed up; however, the biggest gain is 
achieved by using the GMRES scheme. Interestingly, when the GMRES scheme is applied, the advantages of using a 
tridiagonal operator are relatively modest compared to the diagonal (local) operator, which is a very encouraging result 
because in 2-D or 3-D cases an application of anything but a local operator is rather cumbersome. Ng acceleration 
provides a speed-up comparable to the case of no acceleration, but the GMRES scheme is clearly superior. 

To fully appreciate the speed of convergence for all species and energy groups, we present in Figs. |3 — Elthe number 
of iterations required as a function of neutrino energy. The convergence criterion is max((5J/J) < 10^^. For electron 
neutrinos, the convergence is pretty rapid for all energy groups, while it is somewhat slower for the lowest energies for 
the standard ALI scheme without GMRES. With GMRES, the solver converges for most energies in 5—8 iterations, 
for both diagonal and tridiagonal operators. For 9^ neutrinos, the convergence is again slower for low energies. For 
neutrinos, the convergence is slowest for higher energies. This behavior is directly related to the proportion of 
scattering in the deep layers. As is customary in photon transport we define the parameter e as 

e=^. (297) 

Quantity 1 — e is sometimes called a single-scattering albedo. The lower the e, the higher the contribution of scattering, 
and, hence, the transport is more non-local and, thus, numerically more difficult. 

In Fig. El we display e for the three species (full lines for Ve, dotted lines for Ve, and dashed lines for v^^) and for 
three energies (the thickest line for E = 1 MeV, the intermediate for 10 MeV and thin line for 100 MeV). For electron 
neutrinos, e is relatively large for all energies, and thus the scheme converges well. For the lowest energies for 
neutrinos e becomes smaller, which is refiected in slower convergence. For z/^ neutrinos, the lower energies exhibit the 
largest e in the deep layers (close to unity), and consequently the scheme converges very rapidly. For higher energies, 
e becomes very small, which combined with a sharp drop of the source function toward the center, leads to slower 
convergence. Notice that e exhibits a sharp drop beyond 150 - 200 km; this, however, does not lead to any significant 
deterioration of the iteration scheme because these layers are already optically thin. 

Next, we compare the results using three different formal solvers, namely the Discontinuous Finite Element (DFE) 
scheme, a first-order short characteristics (SC) scheme, and the Feautrier scheme. In Fig. [7| the mean intensity J 
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is plotted as a function of radius for four selected energy groups, E = 4.6 MeV, 11.7 MeV, 40 MeV, and 74 MeV. 
Full lines display the DFE results, dotted lines the SC results, and dashed lines the results obtained by the Feautrier 
scheme. Since J spans many orders of magnitude, the lines are almost always indistinguishable on the plots, with 
the exception of J for j/^ neutrinos at higher energies. We, therefore, plot also the relative difference of the mean 
intensities with those computed by the DFE scheme, namely J(solver)/ J(DFE) — 1. The results are displayed on the 
right panels of Fig. [3by dashed (for Feautrier) or dotted (for SC) lines. 

For Ve neutrinos, all three solvers produce mean intensities which are generally within 1% of each other, with the 
exception of the highest energy, where the difference reaches about 9%. However, it should be realized that this 
difference occurs at layers where the value of the mean intensity is some 10 or more orders of magnitude lower that 
the peak value; in view of this fact the overall accuracy is remarkably high. 

For the 9^ neutrinos the accuracy is also pretty high, although somewhat lower than for electron neutrinos (around 
2%); the accuracy decreases to about 10% for the three lower energies around 10 km, and for the highest energy 
around 40 km, where J is very low anyway. The differences between the solvers is largest in the case of neutrinos, 
particularly for the highest energies. Interestingly, the Feautrier solver produces larger differences with the DFE than 
SC at lower energies, while for the highest energies the SC solver is very inaccurate and the Feautrier solver stays at 
the same level of accuracy as at the lower energies. 

In the following tests and production runs, we adopt the DFE solver as our default solver. This is based partly 
on the results displayed above, and partly on the fact that we will adopt DFE for future 2-D simulations, which is 
essentially the only viable choice for irregular or unstructured grids. 

Next, we examine the accuracy of the moment equation solver. We display in Fig. |H1 a comparison of the mean 
intensity computed by the angle-dependent solver (solid lines), and from the moment equation solver. The dotted 
lines represent the original solver, without the sphericity factors and the dashed lines represent the moment solver 
with sphericity factors. As expected, the sphericity factors improve the agreement considerably. In particular, for 
lowest energies the moment equation solver results in a difference of about 10% or more. Using sphericity factors, the 
agreement is within a few percent (except, again, at the highest energy), which is quite reasonable. 

Finally, we compare in Fig. El the full solution with that setting 6 = 0- i.e., assuming isotropic scattering (solid 
lines), and with that setting w = (dashed lines). The differences between the solutions assuming anisotropic and 
isotropic scattering are rather small at small radii, which is quite understandable in view of the large optical depth. 
Farther from the center, the differences become larger, but are still quite modest (typically around 5% for most energy 
groups). Neglecting velocities has a larger effect, reaching about 15% for the lowest neutrino energies. 

In Fig. ^1 we plot a comparison of the total net heating rates for the models displayed in Fig O This quantity 
is germane to the neutrino mechanism of core-collapse supernova explosions. The net heating rate differs relatively 
little, reaching about 8%. In the so-called "gain" region just behind the shock, unlike in the comoving case (Buras et 
al. 2006), adding the velocity-dependent term in the mixed- frame formalism, for which the radiation field is calculated 
in the laboratory frame, leads to additional heating. This is explained simply as the inclusion of the blue-shift of the 
laboratory-frame radiation (both the monochromatic energy and the specific intensity) into the frame of the matter 
due to inward accretion against the outward laboratory-frame flux. This is not to say that the Buras et al (2006) 
analysis is incorrect; the two approaches should give the same results. Rather, it shows that the velocity corrections 
depend upon the frame in which one calculates the radiation quantities. If done in the comoving frame, the numerous 
velocity corrections result in a slight diminution in the heating rate behind the shock, whereas if done in the laboratory 
frame the velocity correction is slightly positive. This has a bearing on the consequences of including the Doppler term 
in VULCAN/2D simulations and we see that doing so will have the opposite effect to that suggested by Buras et al. 
(2006). 

7.2. Implicit coupling: Protoneutron star cooling and deleptonization for fixed density 

We have made several tests of temporal evolution of the initial structure described in t|7.1l At each timestep, we 
perform an implicit update of T and Y^, as described in ^ Since we do not do any hydro in the present tests, we 
keep density and velocity fixed during the evolution. We start with an initial timestep (typically 1 /is), and the next 
timesteps are set dynamically depending on the value of the maximum relative change of T and Y^. Specifically, the 
timestep At is set to At = Ato('5o/|<5max|)^, where Ato is the previous timestep, Jmax is the maximum relative change 
of T or Ye (whichever is larger), and Sq and p are pre-set parameters. We set Sq = 10~^ and p = 1/2 in the following 
tests. 

As discussed in ti6.2l the implicit update of T and Yg is in principle obtained by Newton- Raphson iterations. We 
found that our scheme converges very fast. When we forced the Newton- Raphson loop to determine new T and Ye 
to perform just one single iteration, we found that this procedure leads to results which are indistinguishable from 
letting the iterations proceed until the convergence criterion is reached (typically, the maximum relative change of ST 
or SYe is less that 10~^). The reason is that the second and subsequent iterations typically lead to very small changes 
in the new T and Y^ (and, moreover, only in one or a few radial zones), so this tiny (or questionable) improvement in 
accuracy does not warrant an accompanied increase in computer time. Consequently, in the following tests we have 
forced the number of total iterations of the implicit update of T and Y^ to be 1. 

As in i l7.1l our standard model considers 16 energy groups for each neutrino species, with logarithmically equidistant 
energies between E = 1 MeV and 300 MeV for v,, neutrinos and 100 MeV for other neutrino species. We test this 
setting by considering two models, calculated in serial mode, with 8 groups per species and 32 groups per species, with 
the same energy limits. The results are displayed in Fig. 1111 where we show the temperature, and Fig. 1121 where we 
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show Yf,^ as a function of radius, for the initial configuration, and for 0.3 seconds and 1 second of evolution. To reach 
1 second, the code required between 3500 to 6000 timesteps, with the larger number of timesteps for lower number 
of energy groups. This was offset by a lower computer time per timestep, so finally it required on a 1.6 GHz AMD 
Linux box roughly 590, 730, and 1700 seconds of computer time to reach 1 second of evolution for 8, 16, and 32 energy 
groups per species, respectively. 

Figiires ITTI and IT^ show that the models with 16 and 32 energy groups per species differ only a little, thus validating 
our choice of 16 groups/species as a default model. However, after 1 s of evolution one already sees some differences 
in the temperature (mostly in the core between 10 and 20 km, where the relative difference reaches some 5%), and 
in Ye J where we found differences of also about 5% in the whole range between 70 and 200 km. On the other hand, 
considering only 8 groups per species leads to larger inaccuracies, reaching 40% in T and 10-15% in Y^. 

An important question is whether one has to perform a full angle-dependent radiation transport solution, and thus 
to update the Eddington factor, in every timestep, or whether one may instead update the Eddington factor each n 
timesteps. In Figs. 1131 and ICT we show the results of such a test. As before, we display the temperature and Ye after 
0.3 seconds and 1 second of evolution, updating the Eddington factor every 10, 100, and 1000 timesteps. The results 
are extremely encouraging: they show that updating the Eddington factor even every 1000 timesteps is quite accurate! 
This is not so crucial in the case of 1-D models, where, thanks to a large efficiency of the transfer solver, the timestep 
with the full transfer solution takes only about 2 to 3 times more computer time than the timestep with only the 
moment equation solver. However, this will be very important in 2 or 3 dimensions, where the full angle-dependent 
solver will take proportionally larger chucks of computer time. 

In the case of the models of Figs. El and ^1 it took 930 seconds, 760 seconds, and 730 seconds, to reach 1 second of 
evolution when an update of the Eddington factor was done every 10, 100, and 1000 timesteps, respectively. However, 
we stress that that the code is not fully optimized for 1-D models since we view those models as mere tests for choosing 
the best setups for future 2-D models. 

8. DISCUSSION AND CONCLUSIONS 

In this paper, we develop the mixed-frame formalism for radiation transport in both one and two dimensions and 
provide various solution methodologies and ID numerical tests. Velocity-dependence to 0{v/c), anisotropic scattering, 
energy-dependent cross sections, energy redistribution due to inelastic scattering, and radiation-matter coupling terms 
are derived and incorporated into the algorithm. The equations in cylindrical, spherical, and planar coordinates are 
provided and the effects of bulk and shear radiation viscosity are automatically embedded into the approach. In 
two dimensions, rotation is consistently included, extending the Eddington factor to an Eddington tensor with five 
components and introducing an azimuthal component to the radiation flux. The zeroth- and first-moment equations 
are derived and their roles in the transport solution and the radiation-matter coupling formalism are fully explored. 
Various solution philosophies, such as Discontinuous Finite Element (DEE), Feautrier, and short-characteristics, and 
various convergence approaches, such as the Accelerated-Lambda-Iteration (ALI), Ng acceleration, and GMRES, are 
also implemented and compared. 

Many radiation- hydrodynamic problems do not require an exquisite treatment of spectral line transport, but a 
good treatment of continuum transport. When such is the case, the mixed-frame approach is clearly superior to 
comoving-frame formalisms for which the velocity-dependence is complicated. The virtues of the mixed-frame per- 
spective are many and include simple velocity dependence with no velocity derivatives, straight characteristics, simple 
physical interpretation, and clear generalization to higher dimensions. Moreover, since multi-dimensional radiation- 
hydrodynamic simulations frequently employ Eulerian grids and hydrodynamics, calculation of the radiation quantities 
in the laboratory (Eulerian) frame would seem to be natural. 

We also stress that one can still use the mixed-frame formalism for treating spectral lines, provided one works 
with high enough frequency resolution. In the past, this was considered an excessive requirement. However, modern 
stationary atmosphere models in 1-D already achieve essentially "full" frequency resolution, e.g., with frequency spacing 
of 0.75 fiducial Dopplcr widths throughout the whole spectrum, that employs some 2 to 3x10^ frequency points (Lanz 
& Hubeny 2003). In this spirit, one can use to advantage the mixed-frame approach as well to for such models and 
for general mass motions. 

The equations and algorithms we have developed can be used for the transport of any radiation, in particular photons 
and neutrinos. In addition, the two-moment closure approach to approximate transport that springs naturally from 
our fully angle-dependent formalism is superior to the incorporation of ad hoc flux limiters into the zeroth-moment 
(energy) radiation equation. Not only is the finite "speed of light" effect automatically included in a two-moment 
closure, but the solution can be made arbitrarily accurate, depending upon how often the correct Eddington tensor is 
updated. 

Radiation transport methods are core tools of the astrophysicist. Time-dependent techniques to address transport 
and the coupling of radiation with matter are of central concern to the theorist trying to explain the dynamical 
phenomena of the Universe. Therefore, along with the ongoing advance of the computational arts, the mixed-frame 
formalism and numerical techniques we have derived and tested in this paper should prove of great value for the 
detailed future investigation of complex astrophysical problems. 

We thank Luc Dcssart, Eli Livne, Jeremiah Murphy, and Todd Thompson for fruitful discussions during the course 
of this work and acknowledge support from both the Scientific Discovery through Advanced Computing (SciDAC) 
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APPENDIX 

MIXED-FRAME TREATMENT OF NEUTRINO INELASTIC SCATTERING 

In order to derive the proper expressions for the neutrino inelastic scattering in the mixed frame we need to consider 
very carefully the transformation of the individual components of the transfer equation. Let us first repeat the 
Boltzmann transport equation for the neutrino occupation probability /: 



i| + (n.V)/^C[/], 



(Al) 



where C = Cnes is the collisional integral (net source term) for inelastic neutrino scattering (the dominant contribution 
is typically that of neutrino-electron scattering). Here, we assume that only the inelastic scattering term is present, 
since we have already treated the other terms in 'i'Z.W 
The relation between the occupation probability and the specific intensity is given by eq. Q, 



where the conversion factor h{i') is given by 
The transport equation (|A1|) is written as 



/ = h{v) f . 



1 dl 

c ot 



(A2) 
(A3) 

(A4) 



So far, the equations have been written without any reference to a coordinate frame. Transfer equations ljA4|l and 
HA1|I have the same form in all frames. Since the collision term is evaluated in the frame in which the matter is at rest 
(comoving frame), we will now write the transfer equation explicitly in the comoving frame, viz. 



1 9/o(i^o,no) 
c dto 



+ (no • Vo) Io{vo, no) = b{va)Co , 



(A5) 



where, as usual, subscript indicates quantities in the comoving frame. The left-hand side of this equation transforms 
as 



1 ^ . V. ^ 

-— + no- Vo 

c oU) 



— ^o(t'o,no) 

^0, 



1 ^ / wn' 



n) , 



so that 



-Hn- V ) /(i^,n) = I — ) fe(iyo)Co(t'o,no 



cdt 



We write this equation as 



where 



+n-V )/(:/, n) = C(i/,n) 



C{v, n) = ( — ) h{vQ) Co{vq, no) , 



(A6) 
(A7) 
(A8) 
(A9) 



We now have to express C (z/, n) through the comoving- frame cross sections and redistribution functions and the 
inertial-frame specific intensities. 
The source term Co is given by Bruenn (1985): 



Co(i^o, no) = [1 - faiva, no)] 
-/o(^'o,no) 



rfc^o^^l^ <) <(^o, no • n;,) 



(AlO) 



where i?™'"'" are the scattering kernels. This equation is easily rewritten using the specific intensity: 



0^(70(1^0, no) 



/o(i^o-no) 



6(^o) 



dv' 

d^o^M'^o, nf,) R'oivQ, J^o, no • nf,) 

^0 



/o(t^o,no) 



dio',^ [6(z.^) - loii^i nj,)] i?r (^0, ^^o: no • K) . 

^0 



(All) 
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Because I^/hi, is invariant, i.e. /(fo, no)/6(t'o) = I{v,n)lb{v), we have 



- 1 



1 - 



/o(t/o,no) 



^0 



1 - 



[h{u) - I{v,n)] , 



and, thus, 
where 



C[v, n) = [h{v) - n)] n) - n) [^-(i., n) - n)] 



r?"'=(i/, n) 



(^.,n) = l(^) i I dLo',^bir.',)Rr{r^o,i^ino-n',) 

C ^ .1 .1 l^n 



dv' 

^0(1^0, no) i^o^l'^o, J^o' "0 • no) , 



dv'n 



and 



I duj',^ Io{,.'„, n',) i?r (i^o, '^'0, no • n'„) . 



In evaluating eqs. HA14|I - (|A16|I . we use the following transformation properties 



dv'n 



and 



W,,n',)du',= r-^\ /(^',n')PT du'=r-^\l{u\n')dco' 



dv' 
V 



(A12) 
(A13) 
(AM) 
(A15) 

(A16) 

(A17) 
(A18) 



where we use the primed quantities {v'q^ Hq, v' ^ n') to refer to incoming ("absorbed") neutrinos, while unprimed quanti- 
ties (I'D, no, V, n) refer to outgoing ("emitted") neutrinos. Finally, we assume that the scattering kernels R are expanded 
in Legendre polynomials, and only the first two terms are retained (Bruenn 1985; Thompson et al. 2003): 

RT' = ('^o, 1^',) + (^0, ^',) no-n',. (A19) 



There is an analogous expression for i?o". Although in the comoving frame one has a simple relation between the 
kernels, i?o" — exp[z/o — i'o)/T]i?g"*, we keep both kernels separate because a transformation of the exponential to the 
inertial frame leads to very cumbersome expressions. 
We now absorb the factor 1/c^, as well as the corresponding Legendre factors 1/2 and 3/2, into the kernels and write 



1 



-< - %^{vo, v'o) + <i>i"(^o, i^^) no • nf, . 



(A20) 



I.e. 



2c2 



2c2 1 ' 



(A21) 



and analogously for In the following, we use only the comoving- frame $0 and $1 so we do not need to use a 

more cumbersome, though more consistent, notation such as (<&i)oj where the first subscript (i = 0, 1) would mean the 
order of the corresponding expansion term, and the second subscript (0) would indicate the comoving-frame quantity. 
In order to evaluate eq. (|A14|) . we have to expand $0 and <I>i, 



mvo,v',)^mv,v') 



and 



9$g 
d In V 



n • w 



9$^ 
d In v' 



d\nv dlnv' ^ 

where superscripts "x" stands for "in" or "out" . Equation (|A14p then reads 



(^,n) = (^) J^j^du^'{l-n'-w)I{v',n') 



^o{v, v') — n - w 



9$o 



^i(v,v ) — n ■ w— n • w- 



d In i 



d In v' 



dlnv d In v 

[n • n' — (1 + n • n') (n • w + n' • w)] 



(A22) 
(A23) 



(A24) 
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where we used eq. (^0)1 for transforming no ■ Uq. We express 

ry-(:.,n)=(^) J '^e{v,v',n), (A25) 
where for e we obtain after some algebra that is completely analogous to that used in deriving eqs. (|22|l - 12t)|) : 



According to eqs. 



and, thus, 



Analogously, 



e{v, v' , n) = J{v') 
(?ln v' 



91ni 



( 2$r 



91nz^ 



9 In i^' 



251 and fA26|l . 



(i/, n) = y ^ e(y, n) = (1 - n • w) y ^ e(zy, v' , n) , 



'4 WW) 



d In i^' 



dv' 



9 In I/' 



n) 



■^0 



9 In v' 



^out 



n(j)Out 



9 In i 



out 

1 



Tijrrwk 2$° 



1 (^2$?"' + 



i91n 1/ 



Finally, the "absorption" coefficient is given by eq. (|Af 5f) . 



te\v, n) 



Ko'=(i'o,no) 



(^) i / rfc.^,^5K) [<i>r(^o,^^) + <i>r(K,^^^)no-n[,)] 



(A26) 



(A27) 



(A28) 



(A29) 



(A30) 



where represents the comoving- frame coefficient. Since they do not depend on specific intensity, the integrations 
can be done in advance and the resulting coefficients tabulated. 
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Fig. 1. — The temperature (T), density (p), and the electron fraction (Ye) as a function of radius, for the fictitious model used in our 
numerical tests. 
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Fig. 2. — Convergence pattern for z/g and f/^ neutrinos for E = 8.6 MeV, and for 5 different setups of the iterative transport solver. See 
text for a discussion. 
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Fig. 3. — The number of iterations to rcacii a converged model, defined by the convergence criterion max(5J/J) < 10 ®, for the individual 
setups of the iterative transfer solver, for all energy groups of the i/e neutrino. 
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Fig. 4. — The same as Fig. |2l but for Ue neutrinos. 
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Fig. 5. — The same as Fig. |3] but for i/^ neutrinos. 
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Fig. 7. — A comparison of the three different formal solvers, namely the Discontinuous Finite Element (DFE) scheme, the first-order 
short characteristics (SC) scheme, and the Feautrier scheme. On the left-hand panels, the mean intensity J is plotted as a function of 
radius for four selected energy groups, with E = 4.6 MeV (black), 11.7 MeV (green), 40 MeV (blue), and 74 MeV (red). Full lines display 
the DFE results, dotted lines the SC results, and dashed lines the results obtained by the Feautrier scheme. The right-hand panels display 
the corresponding relative differences of the mean intensities with those computed by the DFE scheme, namely J(solver)/J(DFE) — 1: 
dashed lines - Feautrier; dotted lines - SC. 
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Fig. 8. — A comparison of the relative differences of the mean intensity computed by the angle-dependent solver and by the moment 
equation solver, J(angle dep. solver)/ J(momcnt solver) — 1. The dashed lines represent the original solver, without the sphericity factors; 
solid lines the moment solver with sphericity factors. The color pattern is analogous to that used in Fig. |7| the only difference is an added 
energy group with E = 1 MeV, represented by the light-blue [turquoise] lines. 
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Fig. 9. — A comparison of the relative differences of mean intensity for the full velocity-dependent, anisotropic scattering solution with 
that setting (5 = (that is, assuming isotropic scattering), but keeping the velocity-dependent terms (solid lines); and with setting u = 
(but keeping 5^0 — dashed lines. The color pattern is the same as in Fig. |7| 
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Fig. 10. — A comparison of the net heating rate ("gain") for the full velocity-dependent, anisotropic scattering solution (full lines) with 
that setting <5 = (that is, assuming isotropic scattering), but keeping the velocity-dependent terms (dashed lines); and with setting w = 
(but keeping 5^0)- dotted lines. Crosses (which lie indistinguishably close to the dotted line) represent the solution with both w = 
and (5 = 0. Note that the sign of the tu-correction to the net gain is different from that found in the comoving-frame formalism (Buras et 
al. 2006). See text for a discussion. 
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Fig. 11. — Temperature as a function of radius, for the initial configuration (full black line), and for 0.3 seconds (blue lines/symbols) and 
1 second (red lines and symbols) of evolution. Solid lines display the evolution where we employ 16 energy groups per neutrino species; 
diamonds display the evolution where we employ 32 energy groups per neutrino species, and crosses display models where we employ 8 
energy groups per neutrino species. 
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Fig. 13. — Temperature as a function of radius, for the initial configuration (full black line), and for 0.3 seconds (blue lines/symbols) 
and 1 second (red lines and symbols) of evolution. Diamonds display the evolution where the Eddington factor (/) was updated every 
10 timesteps; crosses display models where the Eddington factor was updated every 100 timesteps, while stars display models where the 
Eddington factor was updated every 1000 timesteps. Notice that in many instances all the symbols overlap. 
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The same as in Fig. 1131 but for Yg. 



